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Preface 


The  purpose  of  this  study  was  to  perform  Monte  Carlo 
simulations  of  neutral  particle  transport  with  primary  and 
secondary  photon  production  as  it  applies  to  a  point 
isotropic  source  within  a  variable  density  atmosphere. 
Specifically,  the  goal  was  to  increase  the  accuracy  of 
neutral  particle  transport  calculations  within  the  atmosphere 
by  accurately  representing  the  variability  of  air  density 
with  altitude.  During  this  study,  analytic  representations 
for  air  density  and  mass  integral,  which  modelled  the 
variation  of  air  density  with  altitude,  were  developed.  The 
generalized  Monte  Carlo  code  MCNP  required  modification  using 
the  newly  developed  atmospheric  model.  The  modifications 
were  thoroughly  tested  before  actual  real  world  simulation 
studies  were  performed. 

Guidance  and  support  for  this  effort  was  received  from 
many  sources.  First,  I  would  like  to  thank  my  thesis  advisor 
LCDR  Kirk  A,  Mathews.  His  expertise  in  the  field  of  radia¬ 
tion  transport  proved  to  be  invaluable  in  enhancing  the 
success  of  my  project.  Second,  I  would  like  to  thank  the 
members  of  my  thesis  committee;  Capt.  Mike  Sabochix  «.id  Ma j . 
Denis  Beller  for  their  assistance  in  the  preparation  of  the 
final  copy  of  this  thesis.  I  would  also  like  to  thank  Judy 
Br iestmeister  and  Tom  Booth  from  Los  Alamos  National  Labor — 
atory  for  their  assistance  with  MCNP.  Finally,  I  would  like 
to  thank  my  wife,  Tina  and  daughter,  Nina  for  the*.r  unending 
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Abstract 

The  generalized  Monte  Carlo  transport  code,  MCNP  was 
modified  for  purposes  of  two-dimensional  neutron-photon 
transport  modelling  in  a  variable  density  atmosphere.  Calcu¬ 
lations  were  performed  using  cylindrical  <r,z)  geometry  and  a 
point  isotropic  neutron-photon  source  at  a  source  altitude  of 
40  kilometers.  Geometrical  dimensions  of  the  problem  were 
large,  approaching  90,000  square  kilometers.  A  collection  of 
neutron-photon  fluence  results  using  ring  detectors  were 
computed  for  field  points  at  constant  ranges  as  a  function  of 
source  elevation  angle  and  for  field  points  at  constant 
altitudes  as  a  function  of  radius.  These  results  may  be 
useful  tc  others  as  benchmark  results  for  further  studies 
using  discrete  ordinates.  Results  showed  a  decrease  in  run 
time  by  a  factor  of  two  was  possible  using  the  modified 
version  of  MCNP.  The  run  time  comparison  was  based  on  a 
study  performed  by  Culp  et  al.  using  the  unmodified  MCNP 
code.  Only  a  fraction  of  the  spatial  cells  required  in  the 
study  by  Culp  et  al.  were  used.  Results  showed  streaming  of 
the  neutron  and  photon  radiation  at  very  close  ranges. 

Buildup  from  multiple  scattering  contributions  proved  to  be 
the  dominant  effect  at  intermediate  ranges.  Attenuation  from 
successive  downscatter  i~  energy  was  shown  to  be  the  dominant 
effect  at  ra  igos  greater  than  about  60  kilometers.  Results 
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indicated  neutrons  to  be  the  dominant  form  of  radiation  based 
on  the  magnitude  of  the  fluences.  Results  obtained  with  MCNP 
were  compared  to  those  obtained  with  the  computer  code  GMAUG- 
II,  which  uses  the  mass-integral  scaling  approximation.  The 
comparison  showed  SMAUG-II  to  be  fairly  accurate  at  points 
close  to  the  source,  but  severe  errors  were  encountered  at 
ranges  exceeding  about  9  kilometers  for  points  at  or  above 


the  source  altitude 
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Executive  Summary 

Introduction 

Computers  have  become  increasingly  powerful  in  recent 
years,  making  two-dimensional  Monte  Carlo  transport  studies 
more  affordable.  Because  the  density  of  air  varies  approxi¬ 
mately  exponentially  with  altitude,  there  is  a  need  to  incor — 
porate  an  efficient,  sophisticated  method  to  represent  this 
density  variation  as  an  alternative  to  using  a  constant  scale 
height  spatial  mesh.  In  a  constant  scale  height  spatial 
mesh,  large  numbers  of  discrete  cells  are  used  to  represent 
the  variation  of  air  density  with  altitude.  Assuming  expo¬ 
nential  behavior,  a  constant  scale  height  factor  (“e-fold" 
distance)  was  used  to  compute  the  density  within  each  region. 
This  large  prol if eration  of  cells  imposed  limitations  on  the 
size  of  the  geometry;  computer  run  time  was  increased  as  a 
consequence  of  using  a  large  number  of  cells. 

If  a  variable  density  atmospheric  model  can  be  imple¬ 
mented  into  a  production  Monte  Carlo  transport  code,  an 
improvement  in  the  overall  accuracy  over  previously  used 
schemes  can  be  achieved.  An  increase  in  the  geometrical 
dimensions  of  the  problem  and  a  decrease  in  run  time  would  be 
an  added  benefit  of  this  density  model. 

A  four-phase  approach  was  used  to  develop,  test,  and  run 
the  modified  MCivlP  code* 
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1  >  analytic  -functions  were  developed  to  represent  both 
the  density  variation  and  cumulative  mass-integral  as  a 
function  of  altitude  based  on  data  extracted  from  the  U.S. 
Standard  Atmosphere  Tables,  1976; 

2)  the  fitting  functions  were  implemented  into  MCNP  and 
the  necessary  auxiliary  modifications  were  accompl ished ; 

3)  the  accuracy  of  the  MCNP  modifications  were  validated 
by  running  a  series  of  Monte  Carlo  simulations.  Results  were 
compared  to  those  obtained  using  the  unmodified  version  of 
MCNP; 

4)  High  altitude  Monte  Carlo  neutron-photon  simulation 
studies  were  performed  using  the  modified  MCNP  code  for  a 
variety  of  test  cases. 

The  following  assumptions  and  approximations  were  used 
in  performing  the  transport  simulations: 

1 )  the  molecular  constituents  of  the  atmosphere 
consisted  only  of  nitrogen,  oxygen,  and  argon  which  normally 
make  up  99.9  percent  of  the  atmosphere; 

2)  the  fractional  composition  of  the  air  remained 
constant  as  a  function  of  altitude; 

3)  there  was  no  time  dependence  in  the  analysis; 

4)  the  source  was  modelled  as  a  point  isotropic  neutron- 
photon  source; 

5)  the  variation  of  air  density  with  altitude  was 
represented  by  a  series  of  continuous  piecewise  e  nential 


interpolating  curves 


Theoretical  Development 


For  tnis  study,  analytic  ■functions  were  developed  to 
represent  the  approximately  exponential  variation  o-f  air 
density  with  altitude  and  the  cumulative  mass-integral  as  a 
■function  of  altitude.  These  expressions  are  given  by 
equations  (2)  and  <6),  respectively.  The  atmosphere  from 
ground  level  up  to  an  altitude  of  1000  kilometers  was  divided 
up  into  36  scale  height  regions.  The  density  variation 
within  each  scale  height  region  was  represented  using  an 
exponential  function.  Each  exponential  function  used  a 
different  density  scale  height  factor  to  represent  the  change 
in  air  density  over  the  scale  height  region.  Density  scale 
height  factors  for  each  region  were  computed  by  fitting  a 
table  of  reference  density  values  using  an  iterative  solver 
routine  employing  Newton's  method.  Mass-integral  data  were 
computed  by  integrating  a  table  of  reference  density  values 
using  a  Simpson's  1/3  Rule  numerical  integration  routine. 

The  data  was  then  tabulated  as  cumulative  distribution  of 
mass-integral  versus  altitude.  The  mass-integral  scale 
height  factors  for  each  region  were  then  computed  using  the 
same  iterative  solver.  A  total  of  72  scale  height  factors 
were  computed  to  16-digj.t  precision  to  enhance  the  accuracy. 

Equation  (2)  computes  air  density  and  equation  (6) 
computes  cumulative  mass-integra 1  as  a  function  of  geometric 
altitude  given  the  appropriate  scale  height  factors.  Equa¬ 
tions  (3)  and  <7)  represent  the  inverted  forms  of  equations 
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Equation  (3)  computes  an  altitude 


} 

(2)  and  (6) ,  respectively, 
given  an  air  density  and  equation  (7)  computes  an  altitude 
given  a  mass-integ  -al .  These  •four  expressions,  together  with 
the  scale  height  factors,  the  corresponding  base  altitudes, 
the  densities,  and  the  mass-integrals,  completes  the  atmos¬ 
pheric  density  model.  The  36  atmospheric  regions  together 
with  the  parameters  for  density  and  mass-integral  are  tabu¬ 
lated  in  Tables  2  and  3,  respectively. 

As  an  alternative  to  Monte  Carlo  computational  modelling 
of  neutron-photon  transport,  faster  computer  codes  utilizing 
the  mass-integral  scaling  approximation  are  used.  Mass- 
integral  scaling  involves  approximating  the  fluence  or  dose 
in  a  variable  density  system  with  previously  computed  results 
of  fluence  or  dose  in  a  homogeneous  system.  The  results  are 
based  on  discrete  ordinates  numerical  solutions  to  the  Boltz¬ 
mann  transport  equation  for  an  infinite  homogeneous  system. 
These  results  are  then  extended  to  a  variable  density  system 
by  using  the  scaling  laws  given  by  equations  (8)  and  <9>  for 
a  cylindrical  <r,z)  geometry.  Since  air  density  varies  with 
altitude,  z,  the  mass-integral  can  be  found  from  equation 
(10).  This  expression  states  that  the  mass-integral  is  a 
function  of  altitude,  density,  and  the  source  elevation  angle 
("look"  angle).  The  geometry  is  shown  in  figure  3. 

The  computer  code  SMAUG-IX,  based  on  mass-integral 
scaling,  relies  upon  discrete  ordinates  calculations  of 
neutron  and  photon  transport  in  a  homogeneous  medium. 
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The  discrete  ordinates  results  are  based  on  a  set  of  least- 


squares  fits  obtained  from  the  computer  code  ANISN,  a 
one-dimensional  discrete  ordinates  code. 

Results  and  Discussion 

Validation  of  the  modified  MCNP  code  showed  good  agree¬ 
ment  as  compared  to  unmodified  MCNP  results.  Comparison 
results  showed  relative  differences  of  approximately  2 
percent  when  using  a  large  number  of  coarse  spatial  cells  in 
the  unmodified  version  of  MCNP.  Relative  differences  of  1 
percent  or  less  were  achieved  when  fine-mesh  spatial  grids 
with  accurate  average  densities  in  each  scale  height  region 
were  used.  Note  that  the  unmodified  MCNP  code,  using  a 
constant  density  spatial  mesh,  produced  results  that  were 
within  2  percent  of  the  continuously  variable  air  density 
model.  However,  a  very  large  number  of  cells  were  required 
in  the  unmodified  MCNP  code.  The  region  of  comparison  in  this 
case  was  restricted  to  radii  not  exceeding  100  km  because  of 
the  excessive  run  time  required. 

Two  basic  simulation  scenarios  were  used  in  performing 
the  Monte  Carlo  transport  studies.  The  first  involved  com¬ 
puting  neutron  and  photon  fluence  values  using  ring  detectors 
at  the  source  altitude  of  40  kilometers.  The  second  scenario 
involved  calculating  neutron  and  photon  fluence  values  using 
ring  detectors  at  constant  ranges  as  a  function  of  source 
elevation  angle.  Ring  detectors  are  usually  necessary  when 
there  are  few  scattering  events  near  the  point  of  interest. 
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Ring  detectors  are  useful  when  there  is  at  least  one  degree 
of  symmetry  about  a  coordinate  axis,  in  this  case,  the 
z-axis.  The  computer  code  SMAUG- 1 1  was  then  used  to 
generate  results  for  the  same  field  point  locations  for 
comparison  purposes. 

Fluence  data  computed  using  MCNP  behaved  in  an  antici- 
paced  manner,  offering  no  particular  surprises.  The  magni¬ 
tude  of  the  neutron  fluence  was  shown  to  be  roughly  an  order 
of  magnitude  higher  than  the  corresponding  photon  results. 
This  indicated  that  neutrons  were  the  principal  form  of 
radiation  in  this  analysis.  This  also  implied  that  the  mean 
free  path  for  neutron  scat'ter  was  short  compared  to  that  for 
photons . 

Along  the  source  co-altitude  band,  streaming  of  the 
neutron  and  photon  radiation  was  observed,  indicating  reduced 
scatter  contributions  at  very  short  ranges  and  low  mass- 
integral.  For  both  neutron  and  photon  radiation,  buildup  at 
intermediate  ranges  from  multiple  scattering  contributions 
increased  for  decreasing  source  elevation  angle  and 
increasing  range.  The  amount  of  buildup  from  the  photon 
radiation  was  less  pronounced  than  that  from  the  neutron 
radiation.  At  long  ranges,  attenuation  of  the  neutron  radia¬ 
tion,  due  to  successive  downscatter  and  diffusion  away  from 
the  points  of  interest,  began  to  dominate  over  buildup.  For 
the  photon  data,  attenuation  from 

successive  scatter  and  absorption  at  low  photon  energies 
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dominated  buildup  at  longer  ranges.  From  3  kilometers  out  to 
at  least  250  kilometers,  the  primary  photon  radiation  was  the 
principal  contributor  to  the  total  coupled  photon  free-field 
-fluence.  Results  indicated  that  the  secondary  photon  radia¬ 
tion  became  prominent  at  very  long  ranges,  and  converged 
towards  secular  equilibrium  with  the  neutron  radiation.  This 
served  as  evidence  that  the  neutron  radiation  acted  as  a 
source  for  the  secondary  photon  radiation. 

With  respect  to  the  approximations  used  in  this  analy¬ 
sis,  SMAUG-II  predictions  of  neutron-photon  transport  charac¬ 
teristics  were  shown  to  be  accurate  to  within  a  few  percent 
at  points  within  a  few  kilometers  of  the  source.  SMAUG-II 
had  the  greatest  difficulty  in  predicting  the  fluence  in 
regions  of  higher  density  and  over  large  mass  ranges. 

Based  on  comparison  results  with  MCNP,  there  were 
several  deficiencies  noted  with  SMAUG-II.  First,  SMAUG-1I 
could  not  account  for  the  escape  of  particles  from  the  prob¬ 
lem,  especially  out  the  top.  Second,  SMAUG-II  is  inadequate 
for  use  at  altitudes  above  the  original  20  kilometer  limit. 
This  stems  from  the  inability  of  StlAUG-II  to  account  for 
multidimensional  transport  effects  such  as  buildup,  downscat- 
ter,  etc.  inherent  in  a  variable  density  atmosphere.  This 
severely  limits  the  usefulness  of  SMAUG-II  as  an  initial 
transport  diagnostic  tool  or  as  a  tool  for  comparison  against 
other  radiation  transport  methods. 

Increasing  the  geometrical  dimensions  of  the  problem 
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over  previously  used  geometries  in  past  Monte  Carlo  studies 
was  now  possible.  The  increase  in  run  time  as  a  consequence 
of  incorporating  a  variable  density  atmosphere  into  MCNP  was 
more  than  offset  by  the  time  saved  in  eliminating  unneeded 
spatial  cells  to  represent  the  variation  of  air  density  with 
altitude.  An  added  benefit  was  to  decrease  some  of  the  user 
overhead  in  preparing  the  input  files. 

In  terms  of  computational  costs,  using  ring  detectors 
proved  to  be  very  expensive  in  comparison  to  cell  or  surface 
tallies.  However,  ring  detectors  were  necessary  because 
large  cells  were  often  used  in  the  problem  geometry, 
rendering  the  use  of  cell  or  surface  tallies  meaningless. 

In  comparison  to  run  times  obtained  using  an  unmodified 
version  of  MCNP,  run  times  from  the  modified  version  of  MCNP 
proved  to  be  extremely  favorable.  The  number  of  cells 
required  was  substantially  reduced,  and  as  a  consequence,  the 
overall  run  time  was  reduced  and  the  problem  dimensions  could 
be  increased. 

SMAUG-II  and  MCNP  neutron  fluence  results  for  the  first 
scenario  are  shown  graphically  in  figures  S  and  6.  Compari¬ 
son  results  between  SMAUG-II  and  MCNP  are  listed  in  Table  4. 
Photon  fluence  results  are  shown  in  figures  7  and  S.  Compar — 
ison  results  are  given  in  Table  6.  SMAUG-II  and  MCNP  neutron 
and  photon  comparison  results  for  the  second  scenario  are 
shown  as  polar  plots  in  figures  9  and  10,  respectively. 
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Summary  and  Conclusions 


The  above  results  Mere  influenced  by  several  -factors. 
First,  the  mean  free  path  for  neutron  scatter  was  short 
compared  to  that  for  photons,  making  neutrons  the  dominant 
form  of  radiation  in  this  analysis.  Second,  results  showed 
that  secondary  photon  transport  was  directly  coupled  with 
neutron  transport,  and  the  photon  radiation  became  increas¬ 
ingly  important  as  range  increased.  Also,  the  primary  photon 
radiation  was  the  primary  contributor  to  thn  total  photon 
free-field  fluence.  Third,  scattering  was  the  dominant 
transport  mechanism  for  both  neutrons  and  photons. 

Recommendations 

Five  specific  recommendations  for  furtner  study  in  this 
area  are  made.  First,  develop  the  modified  MCNP  code  so  that 
it  is  applicable  to  a  wider  variety  of  density-dependent 
transport  problems.  Second,  investigate  alternate  geometries 
and/or  variance  reduction  schemes  to  increase  the  efficiency 
cf  transport  simulations.  Third,  alternative  a-ray-search 
algorithms  might  improve  the  efficiency  of  the  modified  MCNP 
code.  Fourth,  MCNP  should  be  further  modified  to  allow  the 
use  of  multiple  materials  within  the  atmosphere.  Finally, 
modifications  to  allow  the  user  to  select  either  the  variable 
density  version  or  the  unmodified  version  should  be  incor¬ 
porated  into  a  single  code. 


HIGH  ALTITUDE  NEUTRAL  PARTICLE  TRANSPORT  USING 
THE  MONTE  CARLO  SIMULATION  CODE  MCNP 
WITH  VARIABLE  DENSITY  ATMOSPHERE 


I .  Introduction 


Background 

The  most  widely  used  method  today  -for  representing  the 
variation  of  air  density  in  Monte  Carlo  transport  calcu¬ 
lations  is  to  use  spatial  cells,  where  the  density  within 
each  cell  is  constant.  The  variation  of  air  density  with 
altitude  was  then  represented  using  a  series  of  constant 
density  regions. 

Past  studies  performed  by  Shulstad  CIO  in  1976  and  Culp, 
Monti,  Shoemaker,  Wesley,  and  Mathews  C23  in  1990  have  inves¬ 
tigated  different  methods  for  incorporating  a  variable 
density  atmosphere  in  order  to  perform  atmospheric  transport 
calculations.  In  response  to  a  growing  need  for  faster 
running  computer  codes  for  solving  transport  problems  in  the 
1970's,  Shulstad  developed  a  specialized  multigroup  dif¬ 
fusion  code  on  a  two-dimensional  non-orthogonal  coordinate 
system.  Mass-integral  scaling  was  incorporated  into  this 
coordinate  system  in  order  to  model  the  approximately  expo¬ 
nential  variation  of  air  density  with  altitude.  Shulstad's 
efforts  were  successful,  but  were  limited  in  accuracy  by  the 
assumption  of  isotropic  or  at  most  .'inear  isotropic  scat- 
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tering .  Reasonable  accuracy  was  obtained  only  if  the  mean 
free  path  of  the  particle  was  short  relative  to  the  distance 
over  which  the  density  varied  significantly. 

The  advent  of  powerful  computers  in  the  1980 's  and 
1990's  allowed  for  solving  the  transport  problem  using  a  two- 
dimensional  Monte  Carlo  analysis.  The  study  performed  by 
Culp  et  al.  involved  solving  the  transport  problem  of  a  point 
isotropic  source  at  high  altitude  (40  km)  using  Monte  Carlo 
techniques.  The  study  assumed  an  exponential  atmosphere  and 
approximated  the  variation  of  air  density  with  altitude  by 
using  a  series  of  piecewise-constant  spatial  grids  to  repre¬ 
sent  p(z).  The  equation  used  to  represent  this  density 
variation  is  shown  below. 


p-Po©*fl 
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(i) 


where 

p<zA>  *  air  density  at  lower  region  of  i*1’’  altitude  band 
Cg/cm33 

p<Zi.»i>  -  air  density  at  upper  region  of  i*n  altitude 
band  Cg/cm*! 

*  upper  altitude  of  i%*  region  Ckml 
z L  *  lower  altitude  of  i**  region  CkmD 
The  7  km  scale  height  factor  used  in  equation  (1)  represented 
the  vertical  distance  required  for  the  air  density  to 
decrease  by  a  factor  of  e.  This  method  had  limitations 
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because  in  the  real  case  the  atmospheric  density  cannot  be 
represented  by  using  constant  density  altitude  bands  which 
use  a  single  constant  density  scale  height  -factor.  A  contin¬ 
uously  variable  exponential  -function  within  each  region,  each 
with  a  characteristic  scale  height  -factor,  would  be  more 
accurate  -for  representing  the  variation  of  air  density  as  a 
function  of  altitude. 

The  primary  emphasis  of  the  study  by  Culp  et  al . 
centered  on  comparing  mass-integral  scaling  results  from  the 
computer  code  SMAUG-I I  C33  to  Monte  Carlo  results  from  MCNP 
C 4 D .  Results  showed  that  mass-integral  scaling  overestimated 
the  radiation  at  high  altitudes  and  at  long  ranges.  C2i303 

Problem 

Past  studies  have  shown  that  mass-integral  scaling  is 
only  correct  when  applied  to  a  uniform  homogeneous  medium, 
but  is  an  approximation  in  the  real  case.  Likewise,  Monte 
Carlo  calculations  using  a  constant  density  spatial  grid 
throughout  the  problem  is  also  only  an  approximation  to  the 
real  case. 

Because  computers  have  become  increasingly  powerful, 
multidimensional  computational  modelling  of  neutral  particle 
transport  within  the  atmosphere  has  become  much  more  affor — 
dable.  There  is  a  need  to  incorporate  a  more  sophisticated 
variable  density  atmosphere  as  an  alternative  to  using  a 
simple  constant  scale  height  spatial  mesh  or  a  mass-integral 
scaling  approximation . 
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Once  a  variable  density  capability  has  been  incorporated 
into  a  working  production  Monte  Carlo  code,  an  improvement  in 
overall  accuracy  previously  af •forded  by  using  other  atmos¬ 
pheric  modelling  techniques  in  studying  the  transport  problem 
will  be  possible. 

General  Approach 

High  altitude  neutron-photon  transport  studies  within  a 
variable  density  atmosphere  were  based  on  a  four  phase 
approach  as  given  below. 

1)  Develop  analytic  functions  for  obtaining  both  density 
and  integrated  mass  density  (mass-integral)  values.  There 
are  two  functional  forms  for  the  analytical  expressions,  one 
for  density  and  one  for  mass-integral.  These  expressions  can 
be  derived  by  fitting  these  functions  directly  to  a  suffi¬ 
ciently  dense  table  of  atmospheric  density  values  obtained 
from  the  U.S.  Standard  Atmosphere  Tables,  1976  C5D  and 
forming  a  grid  of  density  and  mass-  integral  scale  height 
factors. 

2)  Implement  the  fitting  functions  into  the  Monte  Carlo 
code  MCNP  and  accomplish  the  necessary  auxiliary  modifi¬ 
cations  in  order  to  use  the  derived  analytic  functions  and 
the  scale  height  factors. 

3)  Validate  the  accuracy  of  the  atmospheric  model  by 
running  simple  one-dimensional  test  cases  or  use  simple  input 
files  which  incorporate  a  fine  spatial  grid  with  accurate 
average  densities  in  each  scale  height  region.  Compare 
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results  to  those  obtained  with  the  variable  density  atmos¬ 
pheric  model . 

4)  Perform  Monte  Carlo  simulation  studies  of  neutron- 
photon  transport  within  a  variable  density  exponential  atmos¬ 
phere  using  a  point  isotropic  source  at  high  altitudes  (40 
km)  and  study  the  effects  over  high  mass  ranges 
(>  100  g/cm*).  Compare  MCNP  results  to  those  obtained  with 
SMAUG-1 1 . 


Scope  of  the  Project 

This  study  involved  the  development  of  analytic  repre¬ 
sentations  for  air  density  and  mass-integral  which  modelled 
the  approximately  exponential  variation  of  air  density  with 
altitude.  The  generalized  Monte  Carlo  code  MCNP  required 
modification  so  that,  at  minimum,  it  was  applicable  to  any 
type  of  radiation  transport  study  within  a  variable  density 
atmosphere,  with  some  exceptions.  These  exceptions  are 
explained  in  section  IV.  Specifically  the  above  methods  were 
used  to  generate  a  collection  of  benchmark  results  for  a 
single  source  altitude  at  selected  field  points  (detectors) 
for  two  scenarios:  1)  constant  range  versus  source  elevation 
angle  and  2)  field  point  altitude  bands  versus  range. 

Results  for  the  regions  of  interest  were  mapped  out  in  a 
multidimensional  fashion  giving  fluence  distribution  values 
for  both  range  and  angle.  Fluence  results  obtained  with  MCNP 
were  then  compared  to  those  obtained  with  SMAUG-II.  These 


results  can  be  used  as  benchmarks  for  other  transport  tech- 
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niques  such  as  discrete  ordinates 


Assumptions 

In  order  to  prevent  unneeded  computational  complexity 
which  would  have  little  effect  on  accuracy,  the  following 
assumptions  were  mades 

1)  The  atmospheric  gas  (molecular)  constituents  con¬ 
sisted  only  of  nitrogen,  oxygen,  and  argon  which  normally 
make  up  99.9  percent  of  the  atmosphere} 

2)  The  fractional  composition  of  the  atmosphere  was  not 
a  function  of  altitude; 

3)  There  is  no  time  dependence  in  the  problem,  only 
fluence  calculations  were  considered; 

A)  The  source  was  modelled  as  a  point  isotropic  neutron- 
photon  source.  The  default  spectra  from  the  5MAUG-II  code 
was  used  to  facilitate  the  comparison  between  MCNP  and 
SMAUG-I I ; 

5)  The  variation  of  air  density  with  altitude  was 
modelled  by  using  continuous  piecewise  exponential  curves  in 
each  scale  height  region. 

Sequence  of  Presentation 

Section  II  presents  a  detailed  discussion  of  how  the 
atmospheric  model  was  derived  and  the  sources  of  data  used. 
Section  III  provides  a  general  description  of  the  computer 
codes  used  in  this  analysis,  MCNP  and  SMAUG-I I.  Section  IV 
addresses  the  MCNP  code  modif ications  and  their  implemen- 
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tation.  Section  V  discusses  the  input  parameters  used  in  the 
MCNP  code,  i.e.,  problem  geometry,  materials  used,  selection 
of  tallies,  variance  reduction  techniques,  and  crass  section 
tables.  Section  VI  presents  the  results  and  discussion  of 
the  Monte  Carlo  simulations  as  it  applies  to  the  modifica¬ 
tions  built  into  MCNP.  Section  VI  also  presents  the  methods 
used  to  validate  the  MCNP  modif ic&wions .  Summary  and  conclu¬ 
sions  are  presented  in  section  VII,  which  are  drawn  directly 
from  results  obtained  in  section  VI.  Section  VIII  makes 


recommendations  for  future  work  in  this  area 


II .  Theoretical  Development  - 


Introduction 

Section  I  explained  some  of  the  basic  concepts  of  a 
variable  density  exponential  atmosphere,  and  the  steps  neces¬ 
sary  to  represent  the  density  variation  of  air  with  altitude. 
Section  I  also  presented  the  basic  scope  and  the  assumptions 
used  in  this  study.  This  section  will  present  the  theo¬ 
retical  development  used  in  representing  the  variable  density 
atmosphere  and  a  theoretical  discussion  of  mass-integral 
scaling.  Topics  discussed  includes  the  source  of  atmospheric 
data,  analytic  forms  of  the  density  and  mass-integral 
equations  and  the  determination  of  scale  height  factors  for 
both  density  and  mass-integral. 

Density  Variation  of  the  Atmosphere 

The  variation  of  air  density  with  altitude  is  approxi¬ 
mately  exponential,  that  is,  density  decreases  exponentially 
with  increasing  altitude  up  to  about  110  kilometers.  This 
behavior  is  shown  in  figure  1.  The  exponential  nature  of  the 
air  density  versus  altitude  is  evident  up  to  an  altitude  of 
approximately  110  km.  Above  this  altitude,  the  data  begins 
to  deviate  somewhat  from  a  true  exponential.  The  data  used 
in  figure  1  was  obtained  from  the  U.S.  Standard  Atmosphere 
Tables,  1976  C5i50-73D. 

In  order  to  solve  transp^  problems  within  the  atmos¬ 
phere,  it  was  necessary  to  represent  this  density  variation 
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as  accurately  as  possible  without  greatly  increasing  the 
computational  work. 


Altitude  (km) 


Figure  i  Density  vs.  Geometric  Altitude  from  the  U.S. 
Standard  Atmosphere  Tables,  1976. 


For  early  development  work  of  the  atmospheric  model,  an 
altitude  limit  of  1000  km  was  imposed  because  of  the  1000  km 
limit  in  the  1976  data  tables  C5iSO-733.  In  performing  the 
actual  Monte  Carlo  calculations,  considerable  computer  time 
can  be  saved  by  cutting  the  effective  "top  of  the  atmosphere" 
down  to  an  altitude  where  there  is  much  less  than  one  mean 
free  path  of  atmosphere  remaining  for  either  neutrons  or 
photons.  Above  this  altitude,  the  probability  for  interac- 
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tion  is  very  small.  It  is  often  necessary  to  tailor  the 
limit  of  the  "sensible"  atmosphere  to  a  specific  problem, 
where  high  altitude  detectors  might  be  used. 

Table  1  lists  the  mean  free  path  versus  altitude  and 
density  for  a  2  MeV  neutron,  the  average  energy  in  the  spec¬ 
trum  used  here.  It  is  clear  that  at  about  70  km,  there  is 
less  than  one  mean  free  path  for  neutrons  remaining  within 
the  defined  atmosphere.  There  is  a  certain  probability  that 
contributions  from  particles  that  scatter  downward  from 
collisions  between  the  specified  altitude  limit  and  the  70  km 
suggested  limit  might  be  missed.  Because  of  this,  signi¬ 
ficant  error  could  result  if  a  Monte  Carlo  analysis  is 
performed  using  an  upper  altitude  of  less  than  70  km. 


Table  1:  Mean  Free  Path  for  a  2  MeV  Neutron  vs.  Altitude 


Altitude  CkmD 

Density  Cg/cm33 

MSSBSSBBBtM 

Sea  Level 

-  1 . 2250E-3 

1.23E+4  /  0. :23 

10.0 

4. 1351E-4 

3.63E+4  /  0.363 

20.0 

8.8910E-5 

1.69E+5  /  1.69 

60.0 

3.0968E-7 

4.85E+7  /  483 

70.0 

8 . 2829E-8 

1.81E+8  /  1810 

80.0 

1 . 8458E-8 

8 . 14E+8  /  8140 

90.0 

3.41 60E-9 

4.40E+9  /  44000 

100.0 

5.6040E-10 

2 . 68E+ 10  /  268000 

150.0 

2.0760E-12 

7.24E+12  /  7.24E+7 

200.0 

2.5410E-13 

5.91E+13  /  3.91E+8 

10 


I 


# 


a 


Model lina  the  Density  Variation 

If  the  variation  of  air  density  is  approximated  as 
piecewise  exponential  and  the  "sensible"  atmosphere  is 
divided  up  into  vertical  regions  for  increased  accuracy,  then 
curve-fits  for  the  density  variation  have  the  following 
mathematical  form  (for  Zx  <  2  <  z *♦») 

p  (r)-p(zi)exs|--^~^  <2> 


where 

p(z)  »  density  at  specified  altitude  tg/cm*3 
p(zx>  =  density  at  base  altitude  of  i*H  region  Cg/cm33 
2  =  specified  altitude  Ccm3 
Zi  *  base  altitude  of  i*H  region  [cm3 
S4  *  density  scale  height  factor  of  i*H  region  Ccm3 
The  parameters  needed  in  equation  <2>  to  calculate  density  up 
to  a  maximum  altitude  of  990  km  are  presented  in  Table  2. 
Inverting  equation  (2)  results  in  the  following  expression 
which  calculates  the  altitude  for  a  given  density  (for  p(Zx) 

>  p<Z>  >  p( Zx-l  >  > 

Density  values  were  taken  directly  from  the  U.S.  Standard 
Atmosphere  Tables,  1976  C5i50-733.  To  obtain  good  accuracy, 
the  atmosphere  was  partitioned  into  36  scale  height  regions. 
TK  Solver  Plus  was  used  to  calculate  a  density  scale  height 


1 1 


■factor  -for  each  region  by  forcing 


P  <*<i*i>>  -P  <*i>  exp C- 


given  z&f  z*f  »  ,  pt,  and  pa.fi  to  determine  Sj..  An  iterative 
solver  algorithm  using  Newton's  method  was  used  for  this 
purpose.  The  36  atmospheric  regions  together  with  the  corre¬ 
sponding  parameters  for  density  are  presented  in  Table  2. 


Table  2s  Atmospheric  Parameters  for  Constant  Density 
Scale  Height  Regions 


Region 

Zi_  Ckm] 

Z„  Ckm] 

pi  Ckg/m33 

Sj.  Ckm] 

1 

0 

1 

1.225 

10.30391726068488 

2 

1 

2 

1.1117 

1 0 . 069270550 1 6300 

3 

2 

3 

1 . 0066 

9.83153560749636 

4 

3 

4 

9.0925E-1 

9.60534141126499 

5 

4 

5 

8. 1935E-1 

9.37232106253662 

5 

10 

7 . 3643E- 1 

8.66352196902309 

S3 

10 

15 

4. 135 IE- 1 

6.64086763389816 

6 

15 

20 

1 .9476E-1 

6.3763846529 

9 

20 

25 

8.8910E-2 

6.27630265119 

10 

25 

30 

4 . 00B4E-2 

6.426047540969 

1  1 

30 

40 

1 .8410E-2 

6.545893995756 

12 

40 

50 

3 . 9957E-3 

7.3601 1605750619 

13 

50 

60 

1 . 0269E-3 

8.34195106975831 

14 

60 

70 

3 . 096BE-4 

7.58287361663704 

15 

70 

80 

8 • 2829E-5 

6.66098102745308 

16 

80 

90 

1 . 8458E-5 

5.92758528126849 

17 

90 

100 

3.41 60E-6 

5 . 53227469727463 

18 

100 

110 

5.6040E-7 

5.70413094786334 

19 

110 

120 

9 . 7080E-7 

6.78176406889826 

20 

120 

130 

2 . 2220E-8 

9.97277941535646 

Table  2i  (Continued) 


Region 

Zl.  C kml 

Zh  C km] 

pi  Ckg/m*] 

Sx  Ckm] 

21 

130 

140 

8.1520E-9 

13.2426226009301 

22 

140 

150 

3.8310E-9 

16.3216567546257 

23 

150 

160 

2.0760E-9 

19.194125389299 

24 

160 

180 

1 . 2330E-9 

23.1339249629764 

25 

180 

200 

5. 1940E-10 

27.9741247049562 

26 

200 

240 

2.5410E-10 

34 . 0828528670665 

27 

240 

280 

7.8580E-11 

41.1254588557962 

28 

280 

320 

2.9710E-1 1 

46.8045755271109 

29 

320 

400 

1.2640E-11 

53.1146385350911 

30 

400 

480 

2.8030E-12 

58 . 9065250294756 

31 

480 

560 

7 . 2080E-1 3 

63.6011115831811 

32 

560 

640 

2.0490E-13 

69.8932556471199 

33 

640 

720 

6.5230E-14 

81 .6274032394408 

34 

720 

800 

2.4480E-14 

104.199494345468 

35 

800 

880 

1 . 1360E-14 

141 .881682961411 

36 

8b0 

990 

6.4640E-15 

198.699242408463 

In  Table  2,  the  column  headed  'Region*  contains  the  region 
number,  i,  the  column  headed  'Zi_'  contain  values  of  ztl  and 
the  column  headed  'ZH'  contain  values  of  zx~i. 

The  density  scale  height  factors,  the  base  altitudes  of 
the  regions,  and  the  air  densities  at  the  base  altitudes  were 
tabulated  as  table  look-up  values  and  loaded  into  arrays. 

See  Section  IV,  "Auxiliary  Modifications",  for  more  infor — 
mation . 

Computing  the  Integrated  Mass  Density 

In  order  to  perform  transport  calculations  in  a  variable 
density  atmosphere,  it  was  necessary  to  create  a  table  of 
integrated  mass  density  (mass-integral)  values  as  a  function 
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of  altitude.  It  was  then  necessary  to  derive  a  mathematical 
formula  which  calculated  a  mass-integral  between  any  two 
points  in  space  lying  on  a  vertical  line.  A  correction  for 
points  which  do  not  lie  on  a  vertical  line  is  made  by 
dividing  the  mass-integral  by  the  vertical ,  or  z-direction 
cosine,  using  the  polar  angle  as  a  reference.  A  Microsoft 
QuickBasic  program,  INTEG_MI .BAS,  was  written  which  calcu¬ 
lated  the  mass-integral  by  integrating  the  table  of  density 
values  using  a  Simpson's  1/3  rule  integration.  The  mass- 
integral  data  was  then  tabulated  as  a  cumulative  distribution 
of  mass-integral  versus  altitude.  The  source  code  for  INTEG- 
_M I . BAS  is  listed  in  Appendix  D.  To  increase  accuracy  and 
avoid  round-off  error,  the  mass-integral  data  was  calculated 
in  double  precision  to  16  significant  digits.  In  hindsight, 
it  would  have  been  better  to  integrate  the  mass-integral  from 
the  maximum  altitude  downward  in  order  to  avoid  round-off 
errors.  At  high  altitudes,  the  difference  between  nearly 
equal  mass-integral  values  requires  extreme  precision  in  the 
scale  height  values,  since  their  difference  is  much  less  than 
1.  Integrating  from  the  maximum  altitude  downward  would 
involve  the  difference  between  two  very  different  numbers 
because  of  the  greater  initial  change  in  mass-integral  with 
decreasing  altitude.  Significant  round-off  error  could  still 
be  avoided  even  at  low  altitudes  because  of  the  much  higher 
density  relative  to  the  higher  altitudes,  and  hence  the 
cumulative  mass-integral  would  continue  to  change  signifi¬ 
cantly.  Also,  since  very  few  particles  undergo  transport 
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down  to  altitudes  below  about  20  kilometers,  very  few  calcu¬ 
lations  would  be  required  in  these  regions  relative  to  higher 
altitude  regions.  Figure  2  shows  a  graphical  representation 
of  cumulative  mass-integral  as  a  function  of  geometric 
altitude. 

The  altitude  regions  defined  for  the  density  data  were 
also  used  for  the  mass-integral  data.  An  expression  for  the 
mass-integral  as  a  function  of  geometric  altitude  was 
obtained  by  integrating  equation  \2>  over  each  altitude 
region,  and  has  the  following  mathematical  form  (for  z* 

<  Z  <  Zi*i ) 


AfJ(z)  -MKzJ  + 


J££41 

Si  J 


dz 


<s> 


where 

lil(z)  =  mass-integral  at  specified  altitude  Cg/cm*  3 
MI  (Zt)  =  mass-integral  at  base  altitude  of  i**-*  region 
Cg/cm*  ] 

p(z*>  *=  density  at  base  altitude  of  i*H  region  Cg/cm3] 
z  «  specified  altitude  Ccm] 

Zi  *  base  altitude  of  i*H  region  Ccm] 

Si,  «  mass-integral  scale  height  factor  of  i**  region 
Ccm] 

Performing  the  integration  yields  the  following  expression 
for  the?  mass-integral  as  a  function  of  geometric  altitude 
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Figure  2  Cumulative  Mass  Integral  as  a  Function  o-f 

Geometric  Altitude  for  a  Polar  angle  of  0  Degrees 


MUz)»MI(zi)+p{zi)Si 


The  parameters  needed  in  equation  <6>  to  calculate  the  mass- 
integral  up  to  a  maximum  altitude  of  990  km  are  presented  in 
Table  3.  Inverting  equation  <6>  results  in  the  following 


expression  which  calculates  the  altitude  given  a  mass- 


integral 


1  1  T  (p  <z‘>  sj) 


<7> 


The  36  atmospheric  regions  together  with  the  corresponding 
parameters  -for  mass-integral  are  listed  in  Table  3.  The 
cumulative  mass-integral  at  an  altitude  o-f  990  km,  which 
represents  the  mass-integral  at  infinity,  was  found  to  be 
103S. 635131402448  g/cm2 .  This  value  defines  the  upper  limit 
of  the  atmospheric  model ,  and  was  computed  using  equation 
(6)  . 

Mass-Integral  Scaling  (MIS) 

Describing  radiation  transport  in  detail  involves  using 
the  Boltzmann  transport  equation.  A  full  physics  analytic 
treatment  of  the  Boltzmann  transport  equation  as  it  applies 
to  radiation  transport  is  rarely  attempted.  Because  the 
density  of  air  varies  with  altitude,  numerical  solution  of 
the  Boltzmann  transport  equation  becomes  prohibitive  in  terms 
of  complexity  and  computational  effort.  As  an  alternative, 
radiation  transport  characteristics  are  studied  through  the 
use  of  Monte  Carlo  simulations.  Monte  Carlo  techniques  can 
provide  excellent  results  by  modelling  the  problem  in  two  or 
three  dimensions.  Often,  Monte  Carlo  simulations  can  provide 
results  only  for  specific  cases  because  of  the  high  computa¬ 
tional  cost  associated  with  this  technique.  This  became  the 
primary  reason  for  developing  faster  computer  programs  which 
used  the  mass-integral  scaling  approximation . 
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Table  3:  Atmospheric  Parameters  -for  Constant  Mass- 
Integral  Scale  Height  Regions 


Region 

Zu  C  km ] 

Zh  C km] 

MI  C kg/m3D 

Si.  tkm“»3 

1 

0 

1 

0 . 0000 

10.3412402237278 

2 

1 

2 

1 16.76350000000 

10.0998192414954 

3 

2 

3 

222.60716666667 

9.8669317261947 

4 

3 

4 

318.33433333333 

9.6450353272150 

5 

4 

5 

404.70453333334 

9.4115591457853 

6 

5 

10 

482.43680000000 

8.8809231258504 

7 

10 

15 

763.99446666667 

6.9286R8443723 

8 

15 

20 

91 1 .27230000000 

6.374429797604 

9 

20 

25 

978.75926666667 

6.261929959628 

10 

25 

30 

1009.3796200000 

6.3973690617414 

11 

30 

40 

1023.2862866667 

6.4755778050254 

12 

40 

50 

1032.6629466667 

7.1036242751129 

13 

50 

60 

1034.8067933333 

8.3821 126459594 

14 

60 

70 

1035.4064800000 

7.7650393215337 

IS 

70 

80 

1035.5806097667 

6.83277351 17248 

16 

80 

90 

1035.6241078667 

6.1260854547166 

17 

90 

100 

1035.6332051467 

5.5727092615073 

18 

100 

1  10 

1035.6347923667 

5.6886484922656 

19 

110 

120 

1035.6350561960 

6.198139513  .538 

20 

120 

130 

1035.6351043807 

9.351096091772 0 

21 

130 

140 

1035.6351180274 

12.6721741503841 

22 

140 

150 

1035.6351236653 

15.7852659297959 

23 

150 

160 

1035.6351265031 

18.7015085627943 

24 

160 

180 

1035.6351281111 

22.1902093305744 

25 

180 

200 

1035.6351297362 

27.1428213479546 

26 

200 

240 

1035.6351304712 

32.5722461 137614 

27 

240 

280 

1035.6351310565 

39.9158435804466 

28 

280 

320 

1035.6351312550 

45.8484861083246 

29 

320 

400 

1035.6351313343 

51 .6265976153724 

30 

400 

480 

1035.6351313857 

57.9144605754359 

31 

480 

560 

1035.6351313979 

55.0126655594858 

32 

560 

640 

1035.6351314009 

68.3636856102718 

) 


Table  3*  (Continued) 


Region 

Zu  Ckm3 

Zm  Ckm3 

HI  Ckg/m*3 

C  km"1 3 

33 

640 

720 

1035.6351314019 

78.3404006307542 

34 

720 

800 

1035.6351314022 

98.1437260557627 

35 

800 

880 

1035.6351314023 

139.178054965649 

36 

880 

990 

1035.6351314024 

188.775256150840 

Hass-integral  scaling  is  a  technique  which  involves 
approximating  the  fluence  or  dose  in  a  variable  density 
system  with  previously  computed  results  of  fluence  or  dose  in 
a  homogeneous  system  C2»23.  Discrete  ordinates  numerical 
solutions  to  the  Boltzmann  transport  equation  provide  the 
computed  results  for  an  infinite  homogeneous  system.  The 
computed  results  for  the  homogeneous  system  are  then  extended 
to  a  variable  density  system,  i.e.,  the  atmosphere.  The 
scaling  law  which  relates  the  mass-integrals  in  medium  A  with 
medium  B,  both  with  uniform  density,  is  given  by  C2x23 

P  J!*A  <B> 

where 

p»ra  *  mass  range  at  r»  in  medium  B 
p ArA  =  mass  range  at  r*  in  medium  A 

From  equation  (8),  it  follows  that  the  4nr*  fluence 
present  at  point  A  is  also  present  at  point  B  with  the 
restriction  that  the  density  is  uniform  and  the  same  iso¬ 
tropic  point  source  exists  in  both  systems.  The  scaling  law 
above  can  be  used  to  approximate  the  masr>  range  in  a  homo- 
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geneous  medium  to  the  mass  range  in  a  variable  density  medium 
by  C2:33 


» 


<9> 


P 


P 


P 


P 


P 


i 


P 


where 

p hTh  =  mass-integral  in  a  homogeneous  medium 
rv  B  range  in  a  variable  density  medium 
pv(r'  )  =  position-dependent  air  density 
MI  =  mass-integral 

Figure  3  shows  the  coordinate  geometry  used  in  this 
analysis.  Because  of  rotational  symmetry,  the  azimuthal 
angle  is  not  needed  and  is  not  shown.  Since  air  density 
varies  only  with  altitude,  z,  the  mass-integral  can  now  be 
defined  by  the  following  expression  C2s33 

m-hz)T0» 


where 

2m  *=  altitude  of  the  source 

z*  «  altitude  of  the  detector  (field  point) 
p(z)  *=  altitude-dependent  air  density 
T  =  source  elevation  angle  (see  figure  3) 

For  small  values  of  f,  there  is  little  change  in  alti¬ 
tude,  and  hence  the  density  can  be  calculated  based  on  a 
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Figure  3  Source-Detec tor  Geometry 


simple  average,  <  z.  +  z„)/2.  Note  that  f  represents  the 
source  elevation  angle,  and  thus  a  correction  for  the  mass- 
integral  along  a  slanted  path  (see  figure  3)  has  been  incor — 
porated  into  equation  (10).  For  Monte  Carlo  calculations, 
the  mass-integral  as  computed  by  equation  (6)  is  for  a 
vertical  path  only.  A  correction  for  a  slanted  path  is  found 
by  diviuing  the  difference  in  the  mass-integrals  at  the  two 
altitudes  by  the  z-direction  cosine,  f  (see  figure  3>. 

SMAUG-II  uses  results  of  4nr*  fluence  versus  mass- 
integral  based  on  a  set  of  least-squares  fits  obtained  from 
the  computer  code  ANISN  (ANlsotropic  SN>,  a  one-dimensional 
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discrete  ordinates  code.  More  details  on  SMAUG-II  are 


provided  in  section  III. 


III.  MCNP  and  SMAUG-II  Code  Descriptions 


MCNP  Code 

MCNP  C4:13  is  a  general-purpose,  continuous-energy, 
generalized-geometry ,  time-dependent,  coupled  neutron-photon 
Monte  Carlo  transport  code.  MCNP  is  written  in  FORTRAN  77 
and  was  originally  developed  by  Los  Alamos  National  Labor — 
atory  (LAND.  It  solves  neutral  particle  transport  problems 
and  can  be  run  in  one  o f  three  modest  neutron  transport  only, 
photon  transport  only,  or  coupled  neutron-photon  transport, 
in  which  photons  are  produced  through  inelastic  interactions 
o-f  neutrons. 

The  user  defines  the  problem  parameters  in  a  standard 
input  file,  INP,  which  is  then  read  by  MCNP.  The  input  file 
contains  information  about  the  geometry,  the  source,  the 
materials,  the  cross  sections,  and  the  tallies  used.  A 
sample  input  file  used  in  this  analysis  is  located  in 
Appendix  B.  A  complete  discussion  of  Monte  Carlo  theory,  use 
of  MCNP,  cross  section  libraries,  and  input  cards  is  provided 
in  the  MCNP  user's  manual  C4  3. 

SMAUB-II  Code 

SMAUG-II  C3j  2i4-51  is  a  computer  program  written  in 
FORTRAN  and  was  developed  by  the  Air  Force  Weapons  Labor — 
atory.  It  computes  the  radiation  effects  from  both  neutron 
and  gamma  radiation  from  a  neutron-photon  point  source  within 
the  atmosphere.  Output  results  include  free-field  neutron 
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and  photon  fluences,  energy  fluences,  mean  energies ,  and 
absorbed  tissue  and  silicon  doses  at  field  point  locations 
specified  by  the  user. 

i 

The  source  information  was  provided  by  the  Air  Force 
Weapons  Laboratory  and  was  furnished  with  the  SMAUG-II 
program.  The  spectra  included  a  37-group  neutron  spectrum 
<in  file  TNS.NEV)  and  a  21-group  photon  spectrum  (in  file 
EXP. GAM).  The  default  spectra  used  in  SMAUG-II  was  also  used 
in  MCNP  to  facilitate  comparisons.  The  neutron  and  photon 
spectra  are  given  in  Appendix  F. 

SMAUG-II  employs  mass-integral  scaling  based  on  discrete 
ordinates  (SN)  calculations  of  neutron  and  photon  transport 
charac teristics  in  an  infinite  homogeneous  system.  SMAUG-II 
uses  a  set  of  1711  empirical  transmission  functions  that  are 
least-squares  fits  to  neutron  and  photon  transport  results 
calculated  with  the  computer  code  ANISN  (ANIsotropic  SN). 
ANISN  is  a  one-dimensional  discrete  ordinates  code  for 
performing  transport  calculations  in  a  homogeneous  medium. 
ANISN  uses  the  DLC-31  multigroup  cross  section  set  provided 
by  Oak  Ridge  National  Laboratory  (ORNL)  for  these  calcu¬ 
lations.  ANISN  also  uses  an  air  density  of  1.11  mg/cm*  and 
an  atmospheric  composition  of  79  percent  nitrogen  and  21 
percent  oxygen. 

Additional  information  concerning  SMAUG-II  can  be  found 
in  AFIT  Technical  Report  AFIT/EN-TR-90-S  C2i4-53.  Details  on 
using  SMAUG-II  can  be  found  in  the  SMAUG-II  user's  guide  C33. 
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) 

A  sample  SMAUG-II  input/output  -file  listing  used  in  this 
analysis  is  provided  in  Appendix  C. 
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IV.  MCNP  Modifications 


Once  the  parameters  for  the  atmospheric  model  were 
determined,  modifications  of  MCNP  began.  This  section  will 
present  information  on  MCNP  program  flow  and  an  overview 
concerning  calculations  used  in  the  user-written  subroutines 
and  functions.  Topics  discussed  includei  Program  flow  of 
Subroutines  HSTORY  and  TRANSM,  where  most  of  the  particle 
transport  takes  place,  reference  density,  calculation  of  two 
main  transport  parameters,  and  the  required  auxiliary  modifi¬ 
cations  to  MCNP.  For  additional  information  on  MCNP  program 
flow,  see  the  MCNP  user's  manual  C4*3B-403.  Additional 
specific  information  concerning  the  MCNP  modifications  is 
given  as  an  addendum  to  this  thesis  C71,  where  copies  of  all 
the  affected  MCNP  subroutines,  usei — written  subroutines  and 
functions  can  be  found. 

Transport  Program  Flow 

Particle  transport  in  MCNP  is  controlled  via  the  main 
transport  overlay  MCRUN.  MCRUN  calls  subroutine  TRNSPT  to 
start  particle  histories  via  subroutine  HSTORY. 

Subroutine  HSTORY.  Subroutine  HSTORY  begins  the  parti¬ 
cle  history  by  starting  a  particle  from  the  source  via 
subroutine  STARTP.  All  source  parameters  such  as  initial 
direction,  weight,  and  energy  are  determined.  Initial  source 
parameters  are  next  checked  against  user — supplied  weight  cut¬ 
offs  and  energy  cut-offs  and  some  summary  information  is 
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incremented.  Subroutine  TALLY  is  called  to  score  tally 
contributions  during  a  particle  history.  If  point  detectors 
or  ring  detectors  are  being  used ,  then  subroutine  TALLYD  is 
called  -from  TALLY  to  score  contributions  to  detectors.  If 
DXTRAN  spheres  are  being  used(  then  DXTRAN  is  called  to  score 
OXTRAN  sphere  contributions.  The  weight  window  variance 
reduction  game  is  played  if  requested,  and  energy  splitting 
or  Russian  roulette  is  performed  if  required. 

Back  in  subroutine  HSTQRY,  particle  transport  continues 
by  calling  TRACK  to  determine  the  intersection  of  the 
particle  trajectory  with  each  bounding  surface  of  the  cell. 
The  distances  to  the  nearest  cell  boundary,  DLS,  time  cut¬ 
off,  DTC,  and  DXTRAN  sphere,  DXL,  are  calculated.  Cross 
sections  for  the  cell  are  determined  via  subroutines  ACETOT 
for  neutrons  and  PHOTOT  for  photons.  A  macroscopic  cross 
section,  QPL,  is  calculated  and  then  modified  by  the  expo¬ 
nential  transform  in  subroutine  EXTRAN,  if  necessary.  The 
distance  to  next  collision,  PMF,  is  determined  by  a  random 
number  sampling  of  a  logarithmic  distribution.  Up  to  this 
point,  no  modifications  were  made  to  subroutine  HSTORY . 

The  track  length  in  the  cell  is  determined  by  taking  the 
minimum  of  the  distance  to  collision,  PMF,  the  distance  to 
time  cut-off,  DTC,  the  distance  to  the  nearest  bounding 
surface,  DLS,  and  the  distance  to  a  DXTRAN  sphere,  DXL.  The 
particle  will  undergo  a  collision  only  if  the  distance  to 
collision,  PMF,  is  less  than  the  distance  to  a  surface  cros- 
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sing,  DLS.  Subroutine  TALLY  is  called  to  increment  any  cell 
tallies  if  necessary  and  some  summary  information  is  incre¬ 
mented.  The  particle's  position  is  also  updated.  The 
energies  and  directions  of  the  particles  exiting  a  collision 
site  are  handled  in  subroutines  ACECAS  and  ACEOS. 

Subroutine  TRANSM.  Subroutine  TRANSM  is  called  from 
TALLYD  to  calculate  the  total  transmission  to  the  detector. 
Prior  to  calculating  the  transmission,  subroutine  DDDET  is 
called  to  calculate  the  distance,  00,  to  the  detector  and 
TRACK  is  called  to  calculate  the  distance  to  the  nearest 
intersecting  surface,  OLS. 

Cross  sections  for  the  cell  are  determined  via  subrou¬ 
tines  ACETOT  for  neutrons  and  PH0T0T  for  photons.  Once  the 
microscopic  cross  section  in  the  cell  has  been  determined,  it 
is  multiplied  by  the  atom  density  in  the  cell  to  obtain  the 
macroscopic  cross  section.  The  track  length  in  the  cell,  0, 
is  determined  by  taking  the  minimum  of  the  distance  to  the 
nearest  bounding  surface,  DLS,  and  the  distance  to  the 
detector,  DD.  The  number  of  mean  free  paths  to  the  detector 
or  nearest  bounding  surface,  AMFP,  is  determined  using  the 
values  of  the  macroscopic  cross  section  in  the  cell,  PLE,  and 
the  minimum  distance  to  the  detector  or  nearest  bounding 
surface,  D.  Then  the  locations  of  the  new  cell  and  particle 
position  are  updated. 

Reference  Density 

Modelling  the  variable  density  atmosphere  involved 


28 


defining  a  reference  air  density.  In  the  unmodified  version, 
MCNP  computes  density-dependent  parameters  usinQ  a  discrete 
cell  density  defined  in  the  input  file.  In  the  modified 
version,  the  variation  of  air  density  is  continuous,  and  a 
constant  density  spatial  mesh  in  the  vertical  direction  is  no 
longer  needed}  cells  are  needed  only  for  splitting/Russian 
roulette.  Therefore,  in  the  modified  version  of  MCNP,  the 
air  density  in  each  cell  is  redefined  to  be  the  mass  density 
at  sea-level,  1.225X10-5*  g/cms,  regardless  of  the  values  read 
from  the  input  file.  All  subsequent  calculations  of  density- 
dependent  parameters  will  now  be  calculated  based  on  a  sea- 
level  homogeneous  atmosphere.  At  this  point,  prior  to 
computing  the  distance  to  collision  or  the  transmission  to 
the  detector,  the  MCNP  calculations  are  intercepted.  The 
aforementioned  parameters  are  then  corrected  for  a  variable 
density  atmosphere. 

Prior  to  running  the  MCRUN  overlay,  the  main  overlay 
IMCN  is  run  to  read  the  input  file  to  obtain  the  necessary 
parameters  used  in  defining  the  problem.  Subroutine  NEXTIT 
is  called  to  process  items  from  the  input  file,  i.e.,  to 
store  such  parameters  as  material  density  and  cell  descrip¬ 
tions.  Subroutine  NEXTIT  was  modified  for  the  purpose  of 
redefining  the  density  in  each  cell  to  be  the  density  at  sea- 
level  . 

Transport  Parameters 

In  Monte  Carlo  transport  calculations,  there  are  two 
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principle  parameters  that  involve  determining  the  optical 
path  length  between  two  points  in  spacet  1)  the  sampled 
distance  to  collision  and  2)  the  number  of  mean  free  paths  to 
a  detector  or  nearest  intersecting  boundary.  Both  parameters 
are  closely  related  and  depend  strongly  on  the  atmospheric 
density.  Details  o-f  the  calculation  o-f  each  o-f  these  two 
important  transport  parameters  is  presented  below. 

Distance  To  Collision.  The  distance  to  collision,  PMF, 
computed  in  HSTORY,  is  based  on  the  microscopic  cross  section 
in  the  cell  and  the  atom  density  in  the  cell,  which  was 
rede-fined  in  NEXTIT  to  be  the  air  density  at  sea-level.  This 
is  given  by 

PLE-  TOTM*DEN{  CELL)  <11> 

where 

PLE  *  macroscopic  cross  section  in  current  cell  Ccm-1D 
TOTM  *  microscopic  cross  section  in  current  cell  for 

either  a  neutron  or  photon  reaction  Cbarns/atoml 
DEN(CELL)  =  atom  density  in  current  cell  Catoms/barn-cml 
Since  the  distance  to  collision  is  a  probability  based  on 
density,  it  is  necessary  to  correct  this  value  for  a  variable 
density  atmospherer  Subroutine  EQDIST,  called  from  HSTORY, 
was  written  to  calculate  the  corrected  distance  to  collision 
based  on  values  input  from  HSTORY,  i.e.,  current  particle 
altitude,  z-direction  cosine,  and  distance  to  collision  based 
on  a  sea-level  homogeneous  atmosphere t  Calculation  of  the 
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required  quantities  used  to  determine  the  corrected  distance 
^  to  collision  are  explained  below.  See  -figure  3  for  the 

general  coordinate  geometry  used  in  this  analysis.  The 
variable  names  referenced  in  the  discussion  below  are  the 

*  variable  names  used  in  the  non-MCNP  (user — supplied)  subrou¬ 
tines  or  functions. 

The  unmodified  version  of  MCNP  computes  distances  in 
)  units  of  distance,  centimeters.  However,  the  variable 

density  atmospheric  model  computes  distances  in  units  of  mass 
range,  g/cm* .  Therefore,  before  correction  for  a  variable 
^  density  atmosphere  can  be  made,  it  was  first  necessary  to 

transform  the  MCNP-computed  distance  to  collision,  PMF  (EDST 
in  subroutine  EODIST),  from  units  of  centimeters  to  an  equiv- 

*  alent  mass  range  in  units  of  g/cm*  for  a  sea-level  homo¬ 
geneous  atmosphere.  This  equivalent  distance  was  then 
corrected  for  a  variable  density  atmosphere  using  the  MCNP 

?  modifications. 

The  mass-integral  along  the  particle  path  in  a  homo¬ 
geneous  sea-level  atmosphere  is  calculated  using  the 

*  following  equation 

MIPDH-1 .225x1  O^EDIST  (12) 

*  where 

MIPDH  *  equivalent  mass  range  along  particle  path  in  a 
sea-level  homogeneous  atmosphere  Cg/cm2 1 
'  EDIST  *  distance  to  collision  in  a  sea-level  homogeneous 
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atmosphere  Ccm3 


This  calculation  is  performed  in  -function  DMI. 

The  cumulative  mass-integral,  MICA,  at  the  current 
particle  altitude,  ZO,  can  be  computed  using  equation  (13) 
below  once  the  scale  height  region  has  been  determined.  For 
(z*  <  zO  <  Zi*t ) 


MICA-MI {zj  +p  (Zi>  S 


J.  1-1 

(13) 

•1  < 

H  ^  JJ 

where 

MICA  «  mass-integral  at  current  particle  altitude 
Cg/cm*  3 

MI<Zi)  =  mass-integral  at  base  altitude  of  i**"*  region 
Cg/cm*  3 

p(z*>  *  density  at  base  altitude  of  i**  region  Cg/cm33 
zO  «  current  particle  altitude  Ccm3 
Zi  =  base  altitude  of  i*n  region  Ccml 
Si.  -  mass-integral  scale  height  factor  of  i**  region 
Ccml 

This  calculation  is  performed  in  function  MASI. 

Though  the  collision  altitude,  ZC,  is  still  unknown,  the 
values  of  MIPDH,  MICA,  and  the  z-direction  cosine,  WO,  can  be 
used  to  determine  the  required  cumulative  mass-integral , 

MICP,  at  the  collision  altitude,  ZC.  MICP  is  computed  based 
on  the  known  mass  range  to  collision,  MIPDH.  This  calcu¬ 
lation  is  performed  by  calling  function  MIZ. 


MICP-MICA +MIPDH*  WO 


(14) 


i 

) 


where 

MI CP  *  required  mass-integral  at  collision  altitude 
Cg/cm*  3 

MICA  *  mass-integral  at  current  particle  altitude 
Cg/cm*  3 

MIPOH  =  mass-integral  along  particle  path  to  collision 
in  a  sea-level  homogeneous  atmosphere  Cg/cm* 3 

WO  =  z-direction  cosine 

At  this  point,  a  check  is  made  to  compare  the  computed 
cumulative  mass-integral  at  the  collision  altitude,  MICP, 
with  the  mass-integral  at  infinity,  MI  INF,  defined  to  be 
1035.63S131402448  g/cm* ,  the  cumulative  mass-integral  at  the 
990  km  altitude  limit.  If  MICP  is  greater  than  MIINF,  then 
Subroutine  EQCIST  is  exited  and  PMF  is  set  to  HUGE,  a  very 
large  number  used  by  MCNP,  effectively  eliminating  any  colli¬ 
sions  along  the  current  particle  track.  If  MICP  is  less  than 
or  equal  to  MIINF,  then  the  collision  altitude,  ZC,  is  deter — 
mined  using  the  following  expression 


zc-Zt-SAi 


(MKzJ-MUzc) 


(p(zi)Si) 


US) 


For  cases  where  the  z-direction  cosine,  WO,  becomes  very 
small  (nearly  co-altitude  relative  to  the  current  particle 
position),  the  change  in  density  is  very  small  along  the 
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particle  path.  Therefore,  the  average  density  between  the 
current  altitude  and  the  collision  altitude  is  used.  In  this 


case,  the  difference  between  the  current  altitude  and 
collision  altitude  is  i  10  meters.  The  air  densities  are 
calculated  by  cabling  function  DNTY. 

Finally,  the  distance  to  collision,  PliF  (EDST  in  Subrou¬ 
tine  EUDIST),  corrected  for  a  variable  density  atmosphere, 
can  be  computed  by  calling  function  EQUIVDST.  Function 
EQUIVDST  accepts  as  inputs  the  collision  altitude,  ZC,  the 
current  particle  altitude,  Z0,  the  z-direction  cosine,  WO, 
the  mass-integral  along  the  particle  direction  in  a  homo¬ 
geneous  sea-level  atmosphere,  MIPDH,  and  the  densities 
<calculated  only  if  WO  is  small).  For  cases  where  WO  is  not 
small,  the  corrected  distance  to  collision  is  given  by 

EDST-  J2222L  ci6> 

wo 


where 

EDST  ®  distance  to  collision  corrected  for  a  variable 
density  atmosphere  Ccml 

For  cases  where  WO  is  small,  the  density  is  nearly  constant 
along  the  particle  path,  hence  the  corrected  distance  to 
collision  is  given  by 

EDST- - MIPDH - —  (17) 

[(DENCA+DENCP)  /2] 


•f  V 


where 

DENCA  «=  air  density  at  current  particle  altitude  Cg/cm33 
DENCP  *  air  density  at  collision  altitude  Cg/cm*3 
Mean  Free  Paths  To  Detector  or  Boundary.  Contributions 
to  a  detector  are  made  at  every  source  or  collision  point  by 
creating  and  transporting  a  pseudoparticle  directly  to  the 
detector.  The  total  transmission  to  the  detector  depends  on 
several  factors:  1)  the  exponential  attenuation  through  the 
medium,  2)  a  probability  density  function  for  scatter  toward 
the  detector,  3)  spherical  divergence,  which  accounts  for  the 
solid  angle  effect,  and  4)  the  particle  weight. 

The  only  one  of  the  af o remen tioned  parameters  that 
depends  on  atmospheric  density  is  the  exponential  attenuation 
term.  The  attenuation  term  represents  the  total  attenuation 
of  the  radiation  by  the  atmosphere  along  the  particle  path 
over  a  given  number  of  mean  free  paths.  The  attenuation 
along  the  particle  path  to  a  detector  is  given  by  the 
following  relation 

A-exp(-AMFP)  <1B> 


where 

A  =  attenuation  aiong  the  particle  path 
AliFP  *  number  of  mean  free  paths  to  the  detector  or 
nearest  intersecting  boundary 
At  this  point,  MCNP  computes  the  distance  to  the 
detector,  DD,  using  subroutine  DDDET.  MCNP  next  determines 
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the  minimum  of  the  distance*  to  the  detector  or  the  nearest 


intersecting  boundary.  MCNP  now  computes  the  macroscopic 
cross  section  in  the  cell,  PLE,  for  a  sea-level  homogeneous 
atmosphere  using  equation  (11).  It  is  necessary  to  correct 
the  macroscopic  cross  section  in  the  cell,  PLE,  for  a 
variable  density  atmosphere  using  subroutine  EQDIST.  A  flag 
is  set  prior  to  calling  EQDIST  in  subroutine  TRANSM  to  indi¬ 
cate  the  origin  of  the  calling  statement,  either  subroutine 
HSTORY  or  subroutine  TRANSM.  The  value  of  the  flag,  either  0 
or  1,  will  determine  the  specific  calculations  that  are 
perf ormed . 

Prior  to  calling  EQDIST,  the  reciprocal  of  the  macro¬ 
scopic  cross  section,  XMFP  *  1/PLE,  is  computed  to  give  the 
mean  free  path  in  a  sea-level  homogeneous  atmosphere  (recall 
that  MCNP  first  calculates  all  density-dependent  parameters 
based  on  sea-level  density  in  each  cell).  EQDIST  is  then 
called  from  subroutine  TRANSM  with  the  mean  free  path,  XMFP, 
the  current  altitude,  ZZZ,  the  z-direction  cosine,  WWW,  the 
calculation  flag,  NINP,  and  the  track  length  in  the  cell,  D, 
as  inputs. 

In  subroutine  EQDIST,  it  is  first  necessary  to  transform 
the  MCNP-computed  minimum  distance  to  the  detector  or  nearest 
intersecting  boundary,  D,  from  units  of  centimeter*  to  an 
equivalent  ma*3  range  in  units  of  g/cm?  in  a  sea-level  homo¬ 
geneous  atmosphere.  This  equivalent  mass  range  is  then  used 
to  correct  the  macroscopic  cross  section,  PLE,  for  a  variable 
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density  atmosphere. 

In  subroutine  EQDIST ,  the  altitude  of  the  detector  or 
nearest  intersecting  bounding  surface  is  determined  using  the 
following  equation 

ZD-ZO+WO*DTD 


where 

ZD  *  altitude  of  detector  or  intersecting  boundary  tern! 

ZO  =  current  particle  altitude  Ecml 

WO  =  z-direction  cosine 

DTD  *  distance  to  detector  or  intersecting  boundary  Ccml 

The  mass-integral  along  the  particle  path  in  a  homo¬ 
geneous  atmosphere,  MIPDH,  is  next  calculated  using  equation 
<  12  >  .  Function  DMI  performs  this  calculation. 

Next,  the  cumulative  mass-integral  at  the  current 
particle  altitude,  MICA,  and  the  cumulative  mass-integral  at 
the  detector  or  intersecting  boundary  altitude,  MIDA,  is 
computed  using  equation  (6).  This  calculation  is  performed 
in  function  MASI . 

The  mass-integral  along  the  particle  path  in  a  variable 
density  atmosphere,  MIPD,  can  now  be  determined  using  the 
following  relation  (for  WO  not  too  small) 

JfJPD-  (MIDA* MICA) 

WO 


(20) 


where 

MI PD  =  mass-integral  along  the  particle  path  in  a 
variable  density  atmosphere  Eg/cm*  3 
For  cases  where  the  z -direction  cosine t  WO,  becomes  very 
small  (nearly  co-altitude  relative  to  the  current  particle 
position),  the  average  density  between  the  current  altitude 
and  the  collision  altitude  is  used.  In  this  case,  the  dif- 
-ference  between  the  current  altitude  and  collision  altitude 
is  £  10  meters.  The  air  densities  are  calculated  by  calling 
function  DNTY.  For  small  values  of  the  i-direction  cosine, 
WO,  MI PD  is  computed  using  the  following  relation 

MIPD-DTD*^  j  (21 ) 


where 

DTD  «  minimum  distance  to  the  detector  or  nearest  inter — 
sec  ting  boundary  Ccm3 

DENCA  =  air  density  at  current  particle  altitude  tg/cm*3 
DENCP  *  air  density  at  collision  altitude  £g/cms3 
Recall  that  the  minimum  of  the  distances  to  the  detector 
and  nearest  intersecting  boundary,  D,  was  first  transformed 
into  an  equivalent  mass  range,  M1PDH.  The  corrected  mass 
range  in  a  variable  density  atmosphere,  MIPD,  was  then 
computed.  It  is  now  possible  to  correct  the  mean  free  path 
to  the  detector  or  nearest  intersecting  boundary,  EDST,  for  a 
variable  density  atmosphere.  This  can  be  calculated  using 
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(22) 


\ 

; 


EDST- 


(EDST*MIPDHl 
MI  PD 


where 

EDST  CLHS3  *  corrected  mean  -free  path  Ccm3 

EDST  CRHS3  *=  uncorrected  mean  free  path  Ecm3 

MIPDH  =  mass  integral  along  particle  path  in  a  sea-level 
homogeneous  atmosphere  Cg/cm?  3 

MIPD  «  mass  integral  along  particle  path  in  a  variable 
density  atmosphere  Cg/cm*  3 

Back  in  TRANSM,  the  corrected  macroscopic  cross  section 
in  the  cell,  PLE,  is  found  by  taking  the  reciprocal  of  the 
corrected  mean  free  path,  EDST  (defined  as  XMFP  in  MCNP), 
where  EDST  is  given  by  equation  (22).  Finally,  the  cumu¬ 
lative  total  number  of  mean  free  paths  to  the  detector  or 
nearest  boundary  over  all  cells  traversed  by  the  particle 
trajectory,  is  found  from  the  following  equation 

AMFP-AMFP+  PLE*  D  (23> 

where 

AMFP  CLH53  “  cumulative  total  number  of  mean  free  paths 
to  the  detector  or  nearest  boundary  along 
the  particle  path 

AMFP  CRHS3  *  number  of  mean  free  paths  along  the 
particle  path  in  the  current  cell 
AMFP  here  is  equal  to  PLE*D. 
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Auxiliary  Modifications 

The  auxilliary  changes  necessary  to  complete  the  MCNP 
modifications  are  explained  in  this  section.  Details  of  each 
of  the  affected  subroutines  not  mentioned  above  will  be 
provided  to  complete  the  modification  package. 

Data  Block  AT1NIT.  A  user-defined  data  block  which 
initializes  the  arrays  containing  parameters  which  define  the 
variable  density  atmospheric  model.  The  arrays  are  stored  in 
a  user-defined  common  block  /ATCOM/.  Common  block  /ATCOM/  is 
defined  in  the  file  CM.inc. 

Subroutine  IMCN.  As  mentioned  above,  this  subroutine  is 
the  main  overlay  subroutine,  and  controls  the  input  and 
sorting  of  the  problem  data. 

Modifications:  1)  Comment  printed  to  the  terminal 
informing  the  user  that  the  variable  density  version  of  MCNP 
was  being  run. 

Comdeck  CM.inc.  This  include  file  serves  as  the  main 
common  block  declaration  file  for  all  overlays  used  in  MCNP, 

Modifications:  1)  Common  block  /ATCOM/  defined  and 
contains  the  necessary  arrays  that  store  the  various  para¬ 
meters  defining  the  atmospheric  density  model. 

Comdeck  ZC.inc.  This  include  file  serves  as  an 
auxiliary  comdeck  file  to  comdeck  CM  and  defines  processor — 
dependent  constants  and  parameters  and  also  dimensions 
general  constants  and  I/O  unit  numbers. 

Modifications:  1)  Two  message  lines  defined  to  inform 
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user  of  MCNP  version  and  for  error  statement  in  subroutine 
EQDIST.  2)  Dimension  variables  defining  the  number  of  atmos¬ 
pheric  regions  and  number  of  scale  height  factors. 

Subroutine  NEXTIT.  This  subroutine  controls  the 
processing  of  input  items. 

Modif icationss  1)  Density  for  all  cells  redefined  to  be 
1.225x10*"*  g/cm*,  the  mass  density  of  air  at  sea-level , 
regardless  of  input  value. 

There  are  absolutely  no  changes  required  for  the 
standard  input  file,  INP.  Cell  density  values  defined  in  the 
input  file  are  required,  but  specific  values  are  irrelevant, 
as  they  will  be  redefined  within  MCNP,  and  thus  have  no 
effect  on  results. 
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V.  MCNP  INPUT  PARAMETERS 


This  section  describes  the  various  input  parameters  used 
in  defining  the  problem.  The  discussion  includes  practical 
details  on  the  geometry,  the  source  spectra,  the  material 
specifications,  the  variance  reduction  techniques,  the  cross 
sections,  and  the  tallies  used.  The  above  parameters  are 
defined  in  the  MCNP  generic  input  file,  INP.  A  sample  input 
file  used  in  this  analysis  is  given  in  Appendix  B.  The 
neutron  and  photon  source  spectra  are  given  in  Appendix  F. 

Problem  Geometry 

The  primary  goal  in  determining  the  geometry  for  a 
particular  problem  is  to  model,  as  accurately  as  possible, 
the  geometric  features  and  the  material  properties  that  are 
inherent  in  the  problem  definition.  A  secondary  goal  is  to 
obtain  the  simplest  geometry  possible  without  overly  compro¬ 
mising  its  validity.  By  retaining  certain  symmetries  in  the 
problem,  the  complexity  and  computational  effort  can  be 
reduced . 

The  only  geometric  feature  inherent  in  the  atmosphere  is 
the  airspace  itself.  The  inherent  material  feature  in  the 
atmosphere  is  the  variation  of  air  density  with  changing 
altitude.  Density  values  vary  several  orders  of  magnitude 
over  a  few  tens  of  kilometers,  which  can  drastically  affect 
the  transport  charac teristics  of  the  radiation.  The  goal, 
therefore,  is  to  develop  a  suitable  geometry  which  "contains" 
the  surrounding  airspace  within  certain  geometrical  limits 


while  utilizing  the  variable  density  atmospheric  model. 

Cylindrical  (r,z)  geometry  is  appropriate  because  of  the 
variation  of  air  density  in  only  the  vertical  (or  z) 
direction  and  the  inherent  azimuthal  or  "rotational"  symmetry 
about  the  z-axis.  To  take  advantage  of  this  symmetry, 
cylinders  of  various  radii  were  used  to  define  the  vertical 
surfaces,  while  planes  at  various  altitudes  were  used  to 
define  the  horizontal  surfaces.  A  graphical  representation 
of  the  geometry  used  in  determining  the  neutron  and  photon 
fluences  at  ranges  at  the  source  altitude  of  40  kilometers  is 
shown  in  figure  4. 

As  seen  in  figure  4,  only  the  right  half  of  the  cylinder 
is  shown,  due  to  the  z-axis  symmetry.  The  problem  has  been 
modelled  in  three  dimensions  with  the  restriction  that  rota¬ 
tional  symmetry  be  employed.  The  geometry  has  been  divided 
up  into  discrete  partitions  or  cells.  Each  cell  is  defined 
separately  in  the  input  file  and  referenced  by  a  user — defined 
cell  number.  Each  cell  is  defined  by  referring  to  surfaces 
which  bound  the  cell  by  a  surface  number.  In  this  case, 
surface  numbering  begins  with  the  innermost  cylinder  (surface 
1 > ,  progresses  outward  until  the  outermost  cylinder  is 
reached  (surface  23),  and  then  begins  at  ground  level 
(surface  24)  and  continues  upward  until  the  top-most  surface 
is  reached  (surface  27).  Detailed  information  on  the  basic 
concepts  of  MCNP  geometry  can  be  found  in  the  MCNP  usr-'s 
manual,  Chapter  1,  section  III  C4i 12-193. 
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ALTITUDE,  Z  (km)  SURFACE 


Figure  4  Spatial  Cell  Geometry  for  Co-Altitude  Detectors 


The  cell  structure  shown  in  figure  4  was  chosen  to  provide 
adequate  particle  splitting/Russian  roulette  (see  variance 
reduction  below)  along  source  co-altitude  ranges  out  to  a 
maximum  detector  range  of  250  km.  Variance  reduction  is  not 
needed  outside  this  co-altitude  region.  Since  the  Monte 
Carlo  code  MCNP  was  modified  to  provide  a  continuously 
varying  atmospheric  density  with  altitude,  constant  density 
spatial  cells  are  no  longer  needed  in  the  vertical  direction 
to  represent  the  variation  of  air  density.  Thus,  the  atmos¬ 
pheric  regions  outside  the  co-altitude  region  can  be  repre¬ 
sented  with  just  two  cells.  This  extra  airspace  is  needed 
for  those  particles  that  possess  a  finite  probability  of 
contributing  to  the  final  result  by  scattering  either  upward 
or  downward  and  then  reflecting  towards  the  point  of  inter — 
est.  There  are  only  25  cells  needed  in  the  co-altitude  case 
(see  figure  4)  as  compared  to  181  cells  required  for  the  same 
type  of  study  in  previous  work.  This  considerably  simplifies 
the  problem  geometry  and  saves  computational  effort.  Cell 
geometry  was  appropriately  modified  for  other  source-detec tor 
scenarios. 

There  were  two  source  spectra  provided  with  the  SMAUG-II 
code  C43.  The  spectra  included  37  neutron  energy  groups 
ranging  from  thermal  energies  to  a  maximum  of  14.2  MeV  and  21 
photon  groups  ranging  from  10-3  MeV  to  a  maximum  of  14  MeV. 
The  neutron  and  photon  spectrum  energies  and  corresponding 
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energy  bin  probabilities  are  listed  in  Appendix  F.  The 
default  spectra  provided  with  the  SMAUG-II  code  were  used  in 
MCNP  to  facilitate  the  comparisons  between  SMAUG-II  and  MCNP. 
The  source  spectra  were  not  classified,  but  rather  generic, 
unclassified  spectra  provided  with  the  computer  code 
SMAUG-1 I . 

During  the  course  of  the  Monte  Carlo  simulations,  an 
attempt  was  made  to  provide  for  both  source  neutrons  and 
source  photons  in  the  same  MCNP  run.  As  a  consequence, 
secondary  photon  results  could  have  been  obtained.  To  accom¬ 
plish  this,  it  would  have  been  necessary  to  generate  source 
neutrons  (via  a  neutron  source  spectrum)  and  primary  source 
photons  (via  a  photon  source  spectrum)  within  the  same  MCNP 
run.  These  attempts  were  unsuccessful,  and  it  is  unknown  if 
generating  multiple  particles  from  multiple  source  spectra  is 
possible  within  MCNP.  As  a  result,  each  MCNP  simulation  was 
repeated  twice,  once  using  the  neutron  spectrum  and  once 
using  the  photon  spectrum.  To  obtain  the  total  photon  free- 
field  fluence  at  a  point  of  interest,  the  secondary  photon 
fluence  obtained  based  on  the  neutron  source  spectrum  and  the 
primary  photon  fluence  obtained  based  on  the  photon  source 
spectrum  would  be  summed. 

The  source  spectra  energy  bin  probabilities  were 
provided  with  SMAUG-II,  and  were  also  used  in  the  MCNP  input 
file.  The  energy  bin  probabi 1 i ties  were  normalized  to  one. 

See  Appendix  B  for  an  example  of  a  neutron  source  spectrum 
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definition.  The  units  for  the  spectra  ares  “neutrons  per  MeV 
per  source  neutron"  for  the  neutron  spectrum  and  "photons  per 


MeV  per  source  photon"  for  the  photon  spectrum. 

Material  Specifications 

The  primary  constituents  of  the  atmosphere  are  nitrogen, 
oxygen,  and  argon,  in  order  of  abundance.  These  elements 
normally  make  up  about  99.9  percent  of  the  earth's  atmos- 
phere.  The  remaining  trace  constituents,  which  make  up  the 
remaining  0.1  percent  of  the  atmosphere,  represent  negligible 
neutron  and  photon  absorption  and  scatter  cross  sections,  and 
thus  were  not  included  in  this  analysis.  Specific  material 
specifications  can  be  found  in  the  sample  input  file. 

Appendix  B. 

Variance  Reduction 

The  flexibility  of  MCNP  enables  the  code  to  provide  a 
wide  variety  of  variance  reduction  techniques  to  aid  the  user 
in  obtaining  peak  efficiency  for  a  particular  problem.  In 
addition  to  the  information  provided  in  MCNP  C4 *  122-1483 ,  the 
Los  Alamos  report  "A  Sample  Problem  for  Variance  Reduction  in 
MCNP"  C63  provides  useful  information  concerning  the  appli¬ 
cation  of  variance  reduction  schemes  offered  by  MCNP.  The 
only  variance  reduction  technique  employed  in  this  analysis 
was  geometric  splitting/Russian  roulette,  and  is  discussed 
briefly  below. 

Geometric  SpI ittina/Russian  Roulette.  This  is  the  most 
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widely  used  variance  reduction  technique  in  Monte  Carlo 
analysis.  This  technique  is  the  most  effective  first  step  in 
obtaining  substantial  gains  in  efficiency.  Geometric  split¬ 
ting/Russian  roulette  was  employed  in  this  analysis  to  aid  in 
particle  penetration  over  long  distances  (large  mass  ranges). 
At  low  altitudes,  (ground  level  to  several  tens  of  kilo¬ 
meters),  the  air  density  is  such  that  the  mean  free  path  of 
neutrons  is  on  the  order  of  several  hundred  meters.  If  the 
point  of  interest  is  several  tens  of  kilometers  distant,  then 
this  represents  an  optically  thick  problem,  and  a  large 
amount  of  particle  splitting  is  needed.  Less  splitting  is 
needed  at  the  higher  altitudes  because  of  decreasing  air 
density. 

In  order  to  transport  these  particles  over  long 
distances,  splitting  surfaces  (defined  in  the  problem 
geometry)  are  required.  Geometric  splitting/Russian  roulette 
is  only  needed  in  the  general  direction  of  the  detector 
locations.  The  ultimate  goal  is  to  uniformly  sample  the 
track-length  estimate  of  the  cell  flux  by  maintaining  an  even 
particle  population  in  each  cell  along  the  direction  to  the 
detector.  Additional  details  concerning  this  procedure  can 
be  found  in  the  Los  Alamos  report  LA-103&3-MS  C6i73  or  AFIT 
Technical  Report  AFIT/EN-TR-90-5  C2«153. 

Cross-Sec  t i on  s 

The  primary  evaluated  source  of  cross  section  data  for 
MCNP  has  been  the  ENDF/B  system.  Specific  details  regarding 


data  tables  included  with  MCNP  (version  3B)  and  their  -formats 


are  contained  in  the  MCNP  manual  C4i4851. 

Several  classes  o-f  nuclear  data  tables  are  included  with 
MCNP.  The  principle  nuclear  and  atomic  interaction  data 
classes  (there  are  5  total)  include!  1)  continuous-energy ,  2) 
discrete  energy,  and  3)  multigruup.  Details  on  each  of  these 
data  classes  is  provided  in  the  MCNP  manual.  For  this  analy¬ 
sis,  continuous-energy  data  tables  were  used  throughout  the 
analysis  in  order  to  extract  the  most  accurate  cross  section 
information  from  the  given  data.  For  the  exact  table  used, 
see  the  sample  input  file,  under  material  parameters,  in 
Appendix  B.  A  brief  discussion  of  continuous-energy  data 
tables  is  provided  below. 

Continuous-energy  data  tables  contain  cross  sections  for 
each  type  of  reaction,  given  over  an  energy  grid  sufficiently 
dense  that  linear  interpolation  between  the  data  points 
results  in  an  interpolation  error  of  normally  not  more  than 
one  percent.  Angular  distributions,  photon  interaction  data, 
and  secondary  photon  production  data  are  included  with  this 
type  of  data  table.  Use  of  continuous  energy  data  tables 
will  usually  increase  run  time  and  will  require  a  larger 
amount  of  computer  memory  [4:481 .  The  specific  cross  section 
librarv  used  in  the  MCNP  simulations  was  the  BMCCS  library 
set  taken  from  ENDF/B-IV.  Specific  details  for  this  library 
set  can  be  found  in  the  MCNP  user's  manual  C4 >523-5571. 
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Tallies 


MCNP  provides  a  wide  variety  of  tallies  and  detectors, 
six  types  for  neutrons  and  five  types  for  photons.  Each 
tally  or  detector  result  is  normalized  to  one  source 
particle.  Descriptions  for  these  tallies  are  provided  in  the 
MCNP  manual  C4j80-823.  For  this  analysis,  only  ring  detec¬ 
tors  were  used.  For  details  on  ring  detector  specifications, 
see  the  MCNP  user's  manual  C4i219-246D.  A  brief  discussion 
on  ring  detectors  and  the  reasons  for  using  them  is  given 
below . 

A  point  oetector  is  known  as  a  "next  event  estimator" 
because  it  is  a  tally  of  the  flux  at  a  point  if  the  next 
event  is  a  trajectory  without  further  collision  directly  to  a 
point  detector  C4:B81.  Contributions  to  a  detector  are  made 
at  every  source  or  collision  point  by  creating  and  trans¬ 
porting  a  pseudoparticle  directly  to  the  detector. 

A  ring  detector  tally  is  a  point  detector  tally  where 
the  point  detector  tally  is  sampled  at  a  location  along  a 
ring.  A  ring  detector  is  useful  when  there  is  at  least  one 
degree  of  symmetry  about  a  coordinate  axis,  in  this  case,  the 
z-axis.  Ring  detectors  were  the  ideal  choice  for  obtaining 
time-independent  flux  (fluence)  results  for  both  neutrons  and 
photons  because  they  enhanced  the  efficiency  of  results  over 
those  obtained  with  point  detectors. 

Because  the  use  of  point  defectors  or  ring  detectors 
involves  the  creation  and  transport  of  psuedopartic les ,  they 
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are  computationally  expensive.  However,  ring  detectors 
become  useful  and  even  necessary  when  there  are  few  scat¬ 
tering  events  near  the  detector  location.  In  these  cases, 
the  use  of  cell  or  surface  tallies  is  often  impossible. 

Also,  when  cells  become  very  large,  computing  a  cell  average 
tally  becomes  meaningless,  and  thus  ring  detectors  are  ideal. 
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VI .  Results  and  Discussion 


Introduction 

In  section  IV,  the  implementation  of  the  MCNP  modifi¬ 
cations  for  a  variable  density  atmosphere  was  discussed. 

This  section  will  present  topics  that  include  validation  of 
MCNP  modifications  and  neutron-photon  fluence  results 
obtained  with  MCNP  and  SMAUG-II. 

Validation  of  MCNP  Modifications 

Testing  of  the  MCNP  ir  jdif  ications  was  accomplished  in 
two  phases:  1)  each  subrc  tine  or  function  using  one  of  the 
fitting  equations,  either  density  or  mass  integral,  was 
tested  and  2)  then  the  modified  MCNP  code  was  tested  as  art 
integrated  whole.  Each  phase  of  the  validation  study  is 
presented  in  turn  below. 

Phase  1 .  Each  new  function  or  subroutine  was  indivi¬ 
dually  tested  using  a  personal  computer  and  Microsoft  Fortran 
4.1  to  ensure  calculational  accuracy.  Each  scale  height 
factor  was  computed  using  TK  Solver  to  16  digit  precision. 
This  was  essential  for  two  reasons.  First,  because  MC.nlP 
calculates  variables  in  cgs  units,  the  scale  height  factors 
became  quite  large  <units  of  cm"1).  Slight  truncation  in  the 
scale  height  factor  caused  substantial  errors  in  the  computed 
density  or  mass-integral  data.  Second,  the  cumulative  mass- 
integral  data  was  computed  from  ground  level  upward  to 
infinity.  Therefore,  substantia'  errors  resulted  in  the 
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computed  values  of  the  scale  height  factors  because  often  two 
nearly  equal  mass-integral  values,  especially  at  high  alti¬ 
tudes,  were  involved  in  the  calculations. 

Before  any  of  the  subroutines  were  written,  it  was  first 
necessary  to  verify  both  the  accuracy  of  the  scale  height 
factors  computed  with  TK  Solver  and  the  accuracy  of  the 
fitting  functions.  The  primary  goal  in  performing  the 
testing  was  to  ensure  that  the  computed  values  of  air  density 
closely  approximated  the  reference  values  taken  from  the  U.S. 
Standard  Atmosphere  Tables,  1976. 

For  the  density  and  mass-integral  fitting  functions, 
equations  (2),  (3),  (6)  and  (7)  together  with  TK  Solver  were 
used  to  verify  both  the  density  and  mass-integral  scale 
height  factors  for  each  scale  height  region.  The  density 
scale  height  factors  were  verified  by  choosing  reference 
density  values  at  each  scale  height  region  and  inserting  them 
into  equations  (2)  and  (3).  Once  the  reference  density  at 
the  base  altitude  and  the  altitude  of  interest  were  inserted, 
the  computed  value  of  the  scale  height  factor  could  be  veri¬ 
fied  by  using  the  built  in  iterative  solver  supplied  with  TK 
Solver. 

After  the  scale  height  factors  for  density  were  veri¬ 
fied,  it  was  necessary  to  calculate  the  relative  accuracy  of 
the  computed  density  values.  Using  36  scale  height  regions 
resulted  in  relative  differences,  found  by  comparing  calcu¬ 
lated  density  values  to  reference  density  values,  that  were 
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lees  than  0.1  percent  in  all  but  two  or  three  of  the  upper 
altitude  regions.  The  maximum  relative  difference  was  found 
to  be  less  than  1  percent  over  all  altitude  ranges.  In  the 
upper  altitude  regions,  the  size  of  the  scale  height  regions 
increased,  leading  to  larger  and  less  accurate  scale  height 
factors.  These  upper  atmospheric  regions  were  not  subdivided 
into  smaller  regions  in  order  to  decrease  the  computer  over — 
head  associated  with  the  array  searches  and  because  the  upper 
atmospheric  regions  are  far  less  important  than  the  inter¬ 
mediate  altitudes  which  are  near  the  source  altitude. 

The  cumulative  mass-integral  data  was  computed  by  inte¬ 
grating  a  table  of  reference  density  values  using  a  Simpson's 
1/3  Rule  numerical  integration  algorithm  written  by  the 
author  (see  Section  IV,  "Computing  the  Integrated  Mass- 
Density").  Simpson's  1/3  Rule  method  has  a  global  error  of 
0<tv*>,  where  h  is  the  vertical  panel  width  for  the  altitude 
data.  The  errors  associated  with  such  a  high  order  numerical 
integration  routine  are  small,  and  on  the  order  of  the 
relative  errors  associated  with  the  density  results.  The 
mass-integral  scale  height  factors  were  verified  in  much  the 
same  way  as  the  density  scale  height  factors.  In  this  case, 
mass-integral  values  were  chosen  at  random  for  each  scale 
height  region  instead  of  reference  density  values. 

Fast  studies  performed  by  Culp  et  al .  relied  on  a 
constant  density  spatial  mesh  using  a  single  7  kilometer 
scale  height  factor  to  represent  the  variation  of  lir  density 


with  altitude.  The  errors  introduced  in  representing  the 
density  variation  in  this  way  could  be  substantial,  ranging 
■from  just  a  few  percent  -for  narrow  altitude  regions  to  SO 
percent  or  more  for  wide  altitude  regions.  Narrowing  the 
altitude  regions  increased  the  relative  accuracy  of  the 
density  variation,  but  at  a  cost  of  increased  run  time  due  to 
the  increase  in  the  number  of  cells  needed  in  defining  the 
problem. 

Several  iterations  were  required  to  determine  the  proper 
number  of  scale  height  regions  (36)  in  order  to  obtain  the 
lowest  possible  relative  error  (<  1  percent).  The  scale 
height  factors  for  both  density  and  mass-integral  were  recal¬ 
culated  each  time  the  boundaries  of  the  scale  height  regions 
were  changed. 

To  fully  model  the  density  variation  in  the  atmosphere, 
subroutines  and  functions  were  written  which  were  later 
implemented  as  auxiliary  modules  in  the  MCNP  code.  These 
subroutines  and  functions  are  explained  in  Section  IV.  The 
parameters  for  the  variable  density  exponential  atmosphere, 
which  included  the  altitudes,  the  density  values,  the  compu¬ 
ted  mass  integral  values,  and  the  scale  height  factors,  were 
loaded  into  data  arrays  as  table  look-up  values.  Each 
subroutine  or  function  module  was  tested  separately  using  a 
personal  computer  and  Microsoft  Fortran  4.1  to  ensure  calcu- 
lational  accuracy.  Calculations  for  both  the  density  and 
mass-integral  equations  presented  in  Section  IV  were  first 
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performed  by  hand .  Calculations  were  then  verified  by 
computer  before  incorporating  them  into  a  particular  subrou¬ 
tine.  Next,  each  subroutine  was  individually  tested  for 
calculational  accuracy. 

Phase  2.  Once  all  the  subroutines  had  been  verified  for 
proper  operation,  they  were  incorporated  into  the  MCNP  code 
and  then  tested.  A  two-part  approach  was  used  to  validate 
the  accuracy  of  the  modified  MCNP  code  results!  1)  run  a 
Monte  Carlo  simulation  from  a  past  study  (Culp  et  al.)  and  2> 
run  a  series  of  simpler  MCNP  simulations  using  a  scaled-down 
geometry. 

Phase  2.  Part  i.  The  study  chosen  was  an  investi¬ 
gation  of  neutron  fluences  using  ring  detectors  at  surfaces 
of  constant  range  versus  selected  source  elevation  angles 
C2:291.  A  maximum  detector  range  of  SO  kilometers  was 
imposed  because  of  the  number  of  cells  required  when  running 
the.  unmodified  MCNP  code.  Fluence  results  from  the  previous 
study  were  compared  to  those  obtained  with  the  modified  MCNP 
code.  The  comparison  showed  a  relative  difference  of  approx¬ 
imately  2  percent.  As  a  first  step,  this  comparison  demon¬ 
strated  that  the  results  obtained  with  the  modified  MCNP  code 
were  reliable. 

Phase  2.  Part  2.  As  a  next  step  in  the  validation, 
simpler  simulation  studies  using  the  unmodified  MCNP  code 
were  performed.  These  studies  were  accomplished  with  scaled- 
down  geometries.  Fine-mesh  spatial  grids  with  accurate 
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Average  densities  in  each  scale  height  region  were  incor¬ 
porated  into  the  input  files.  The  average  densities  in  each 
region  were  computed  using  the  newly-developed  atmospheric 
density  model  to  give  the  most  accurate  comparison.  The 
simulations  involved  small  geometries  using  cylinders  of  a 
few  tens  of  kilometers  in  radius  and  height.  The  atmospheric 
region  was  centered  around  the  primary  source  altitude  of  40 
kilometers.  Again,  fluence  results  using  both  ring  detectors 
and  cell  tallies  from  the  unmodified  MCNP  code  were  compared 
with  those  obtained  with  the  newly-modified  MCNP  code.  The 
comparison  showed  a  relative  difference  of  approximately  1 
percent.  A  closer  comparison  was  reasonable  since  the 
density  variation  as  modelled  using  a  spatial  grid  more 
closely  approximated  a  continuously  varying  density.  Several 
validation  simulations  of  this  type  were  run,  using  ring 
detectors  and  cell  tallies  at  varying  locations  out  to  a 
maximum  range  of  50  kilometers.  The  maximum  relative  dif¬ 
ference  between  the  unmodified  and  modified  MCNP  codes  was 
found  to  be  not  more  than  1  percent.  This  also  demonstrated 
that  the  modified  MCNP  was  functioning  properly. 

To  increase  confidence  in  the  proper  functioning  of  the 
modified  MCNP  code,  comparisons  were  made  of  the  diagnostic 
information  printed  in  the  MCNP  output  file.  The  diagnostic 
information  was  compared  to  see  if  hidden  bugs  could  be 
found.  Important  parameters  for  comparison  were t  1)  particle 
tracks  escaping,  number — weighted  energy,  f lux-weighted 
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energy,  and  detector  diagnostic  information  such  as  detector 
transmissions.  For  details  on  interpreting  MCNP  output,  see 
the  MCNP  user's  manual  C4* 312-4133. 

MCNP  and  SMAU6-I 1  Results 

MCNP  neutron  and  photon  fluence  results  are  listed  in 
Appendix  A  together  with  corresponding  error  estimates.  The 
SMAUG  neutron  and  photon  results  are  also  listed  in  Appendix 
A.  Only  the  most  representative  MCNP  and  SMAU6-II  data  have 
been  presented  here.  Appendix  A  contains  numerous  other  MCNP 
and  SMAUG-II  results  for  other  al ti tude-range  combinations. 

Neutron  Fluence  at  40  km.  Figure  5  shows  a  comparison 
between  MCNP  and  SMAUG  neutron  fluence  results  at  an  eleva¬ 
tion  of  40  kilometers  (source  altitude)  versus  radius  (slant 
range  in  this  case).  Figure  6  shows  the  same  results  versus 
mass-integral .  The  MCNP  results  were  computed  using  ring 
detectors  and  all  results  are  normalized  to  one  source 
particle.  For  comparison  purposes,  a  line  showing  spherical 
divergence  for  one  source  particle  has  been  plotted. 

Close  to  the  source  (0-3  kilometers),  the  neutrons 
exhibit  a  streaming  (spherical  divergence)  behavior.  At  3 
kilometers,  the  mass  range  is  not  more  than  1.2  g/cm* ,  as 
seen  in  figure  6.  At  intermediate  radii  (3-60  kilometers), 
buildup  from  neutron  scattering  contributions  dominates  over 
attenuation  along  the  source  co-altitude  band.  This  results 
in  a  higher  neutron  fluence  than  that  predicted  by  spherical 
divergence.  At  these  ranges,  the  mass  range  is  between  1.2 
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Co-Altitude  vs.  Radius  (Slant  Range) 


g/cm*  at  3  kilometers  and  24  g/cm*  at  60  kilometers  as  shown 
in  figure  6.  At  greater  ranges  (60-250  kilometers),  attenua¬ 
tion  of  the  neutron  radiation  begins  to  dominate  over  build¬ 
up.  This  results  in  a  lower  neutron  fluence  than  that  pre¬ 
dicted  by  spherical  divergence.  At  these  ranges,  the  mass 
range  is  24  g/cm*  at  60  kilometers  and  99.9  g/cm*  at  250 
kilometers.  The  most  important  transport  mechanisms  which 
account  for  this  behavior  are  explained  below  for  the  MCNP 
results . 

1)  At  short  ranges  <0-3  kilometers),  the  relatively  low 
air  density  at  this  altitude  leads  to  a  small  accumulated 


59 


Figure  6  MCNP  vs.  SMAUG  Neutron  Fluence  at  Source 
Co-Altitude  vs.  Mass-Integral 

mass  range.  Therefore,  neither  significant  attenuation  due 
to  energy  loss  of  the  neutrons  nor  buildup  from  scattering 
contributions  has  yet  occurred. 

2)  At  intermediate  ranges  (3-60  kilometers),  the  mass 
range  has  increased  to  the  point  where  buildup  is  occurring. 
In  this  case,  the  neutrons  are  scattering  and  contributing 
several  times  to  the  free-field  fluence  at  the  point  of 
interest.  Because  the  scatter  cross  sections  for  neutrons 
generally  follow  a  1/fE  <l/v)  relationship,  the  higher  energy 
neutrons  are  scattering  down  in  energy,  and  hence  the  proba- 
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bility  for  an  additional  scatter  increases.  The  loss  of 
particles  due  to  successive  scatters  and  their  subsequent 
loss  of  energy  is  more  than  compensated  for  by  the  scattering 
contributions.  This  is  the  reason  for  the  buildup  at  inter¬ 
mediate  ranges. 

3)  ^t  some  point,  the  loss  of  neutrons  due  to  their 
decreasing  energy  overwhelms  the  contribution  to  buildup  from 
neutron  scattering  events,  and  attenuation  begins  to  aominate 
over  buildup.  This  occurs  at  ranges  between  60  and  250  kilo¬ 
meters,  where  the  free-field  fluence  just  equals  the  value 
predicted  by  spherical  divergence.  The  buildup  factor  at 
this  point  is  unity.  Beyc  d  60  kilometers,  attenuation 
dominates  over  buildup  from  scatter  because  the  neutrons  no 
longer  have  enough  energy  to  further  contribute  to  the 
fluence  via  scattering  contributions. 

The  same  general  trends  seen  in  the  MCNP  results  can  be 
seen  in  the  SMAUG  results  shown  in  figures  5  and  6.  Close  to 
the  source,  MCNP  and  SMAUG  show  good  agreement.  At  ranges 
greater  than  about  9  kilometers,  SMAUG  begins  to  severely 
overpredict  the  neutron  fluence  relative  to  the  MCNP  results. 
Table  4  contains  the  ratios  of  SMAUG  to  MCNP  neutron  fluence 
results  versus  range  along  the  source  co-altitude  band. 

Close  to  the  source,  the  SMAUG  fluence  is  somewhat  higher, 
but  the  ratio  is  relatively  close  to  1,  which  confirms  the 
accuracy  of  the  MCNP  results  at  close  range.  At  9  kilo¬ 
meters,!  SMAUG  begins  to  overestimate  the  fluence  signifi- 
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Table  4t  Ratio  of  SMAUG  to  MCNP  Neutron  Fluence 
at  Source  Co-Altitude  versus  Range 


cantly  relative  to  MCNP,  where  the  ratio  is  1.12.  The  ratios 
generally  increase  with  increasing  range  out  to  a  maximum 
ratio  of  14.11  at  a  range  of  250  kilometers.  The  most  impoi — 
tant  reasons  which  account  for  the  behavior  of  the  SMAUG-II 
results  are  stated  below. 

1)  SMAUG-II  performs  calculations  of  mass-integral 
scaling  based  on  approximate  fits  to  one-dimensional  discrete 
ordinates  solutions  of  neutron  and  gamma  transport.  Mass- 
integral  scaling  is  strictly  correct  only  in  a  uniform  homo- 
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geneoun  medium,  and  works  relatively  well  only  when  the  mean 
-free  path  of  the  particle  is  short  compared  to  the  distance 
over  which  the  density  varies  significantly.  In  this  analy¬ 
sis,  the  atmosphere  does  vary  significantly  in  density  over  a 
relatively  short  vertical  distance,  and  hence  errors  result¬ 
ing  from  the  mass-integral  scaling  approximation  are  expected 
on  this  basis  alone. 

2)  The  mass-integral  scaling  approximation  used  in 
SMAUG-II  assumes  that  geometric  divergence  and  scattering  can 
be  separated.  Therefore,  SMAUG-II  would  predict  that  the 
sa  e  buildup  factor  would  exist  for  two  detectors,  one  above 
and  one  below  the  source,  provided  that  the  slant  ranges  to 
each  detector  would  result  in  both  4nR*  fluences  being  equal. 
A  two-dimensional  transport  calculation  by  MCNP  would  predict 
that  the  radiation  would  undergo  less  buildup  toward  the 
upper  detector  than  the  lower  detector.  This  is  due  to  the 
lower  air  density  above  the  source,  causing  fewer  scatters 
and  more  particle  leakage  out  the  top  of  the  problem,  which 
decreases  the  contribution  to  the  free-field  fluence.  At 
high  altitudes,  the  air  is  tenuous  and  the  probability  for  a 
particle  escaping  the  problem  increases  because  there  are 
fewer  scatters  in  this  region.  At  very  low  altitudes,  the 
mean  free  path  for  scatter  is  short,  and  thus  the  particle 
undergoes  many  scatters  along  its  tortuous  path.  This 
results  in  a  successive  loss  of  particle  energy.  As  a  conse¬ 
quence,  the  -article  diffuses  away  from  the  point  of  interest 
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which  prevents  contribution  to  the  free-field  fluence.  In 
this  case,  SMAUG-II  would  tend  to  overpredict  the  fluence 
everywhere,  with  the  overprediction  becoming  more  severe  as 
the  mass  range  increases.  Therefore,  SMAUG-II  would  not  be 
expected  to  account  for  simple  particle  leakage  from  the 
problem.  The  results  shown  in  figure  5  shows  this  to  be 
true. 

The  mechanism  of  "over  the  top"  contribution  of 
particles  initially  travelling  into  the  upper  elevations, 
undergoing  collision,  and  then  reflecting  downward  and  con¬ 
tributing  to  the  final  result  at  long  ranges,  was  not 
observed  in  these  results.  A  mechanism  of  this  kind  is  not 
accounted  for  by  SMAUG,  and  thus  if  this  mechanism  were 
significant,  SMAUG  would  tend  to  underestimate  the  fluence  at 
these  ranges.  Clearly,  results  shown  in  figure  5  do  not 
reflect  this. 

Overall,  these  results  behaved  in  an  anticipated  manner, 
but  figures  5  and  6  clearly  show  the  deficiencies  in  SMAUG- 
1 1 . 

Photon  Fluence  at  40  km.  Figure  7  shows  a  comparison 

a 

between  MCNP  and  SMAUG  photon  f'uence  results  at  an  elevation 
of  40  kilometers  (source  altitude)  versus  radius  (slant  range 
in  this  case).  Figure  8  shows  the  same  results  versus  mass- 
integral.  The  MCNP  results  were  computed  using  ring 
detectors  and  al *  results  are  normalized  to  one  source 
particle.  For  comparison  purposes,  a  line  showing  spherical 
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Co-Altitude  vs.  Radius  (Slant  Range) 


divergence  for  one  source  particle  has  been  plotted. 

Both  the  primary  and  secondary  photon  fluences,  computed 
from  MCNP ,  have  been  platted  separately  in  order  to  chart  the 
behavior  of  each  and  to  see  in  what  region  each  is  important 
in  contributing  to  the  total  fluence.  In  performing  these 
calculations ,  one  source  (primary)  photon  was  emitted  for 
every  source  neutron  emitted.  This  is  not  necessarily  repre¬ 
sentative  for  every  source;  a  detailed  analysis  to  determine 
the  optimum  ratio  of  source  photons  to  source  neutrons  was 
not  performed  due  to  time  constraints. 

Figure  7  clearly  shows  that  the  general  behavior  of  the 
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Co-Altitude  vs.  Mass-Integral 

photon  radiation  generally  parallels  the  behavior  of  the 
neutron  data  shown  in  figure  5.  At  ranges  close  to  the 
source,  the  primary  photon  radiation  exhibits  a  streaming 
nature.  At  intermediate  ranges,  buildup  from  primary  photon 
scattering  contribu -ions  dominates  over  attenuation  from 
absorption  and  energy  losses  from  scattering.  Buildup  in 
this  case  is  less  pronounced  than  for  neutrons  due  to  the 
difference  in  interaction  mechanisms  and  their  associated 
cross  sections.  At  long  ranges,  attenuation  dominates  over 
buildup  and  the  fluence  is  below  that  predicted  by  spherical 
divergence.  The  contribution  to  the  free-field  photon 


fluence  is  dominated  at  all  ranges  by  the  primary  photon 
radiation.  Figure  7  also  shows  that  the  relative  difference 
between  the  primary  and  secondary  photons  is  not  constant. 

The  relative  difference  in  fact  decreases  with  increasing 
distance  (mass  range).  Table  S  contains  the  ratios  of  SMAUG 
to  MCNP  total  (primary  plus  secondary)  photon  fluence  results 
versus  range  along  the  source  co-altitude  band.  Close  to  the 
source,  the  SMAUG  photon  fluence  agrees  well  with  the  MCNP 
fluence  results.  In  this  region  the  ratio  is  relatively 
close  to  1,  which  confirms  the  accuracy  of  the  MCNP  results 
at  close  range.  At  18  kilometers,  SMAUG  begins  to  over — 
estimate  the  fluence  significantly  relative  to  MCNP,  where 
the  ratio  is  1.11.  The  ratios  generally  increase  with  in¬ 
creasing  range  out  to  a  maximum  ratio  of  6.21  at  a  range  of 
250  kilometers.  Reasons  for  the  behavior  of  the  photon 
radiation  computed  by  MCNP  are  given  below. 

1)  Photons  interact  differently  than  do  neutrons. 

Photons  can  undergo  photoelectric  effect,  with  a  cross 
section  which  varies  approximately  as  1/E3,  Compton  scatter, 
with  a  cross  section  which  varies  little  with  energy,  and 
pair  production,  with  a  cross  section  which  increases  rapidly 
with  increasing  energy.  Intuitively,  differences  in  the 
transport  mechanisms  would  be  expected  based  on  these  inter¬ 
action  mechanisms.  However,  results  in  figure  7  show  that 
the  behavior  of  the  photon  radiation  follows  the  same  general 
pattern  of  streaming,  buildup,  and  attenuation  as  the  neutron 
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Table  St  Ratio  of  SMAUG  to  MCNP  Photon  Fluence 
at  Source  Co-Altitude  versus  Range 


Range  Ckml 

Ratio 

3 

1.03 

6 

1 .01 

9 

1 .01 

12 

1 .03 

15 

1.07 

18 

1.11 

25 

1 .30 

50 

2.25 

75 

3.39 

100 

4.45 

125 

5.39 

150 

5.88 

175 

6.10 

200 

6.60 

225 

6.50 

250 

6.21 

radiation.  The  streaming  behavior  is  predicted  based  on 
similar  reasons  given  for  neutron  data.  Both  the  air  density 
and  the  accumulated  mass  range  are  low,  therefore  neither 
significant  attenuation  nor  buildup  from  photon  scattering 
contributions  has  yet  occurred.  Buildup  occurs  at  interme¬ 
diate  ranges  because  both  range  and  mass  range  have 
increased,  resulting  in  more  photon  scatters  which  leads  to 
multiple  contributions  to  the  free-field  fluence.  However, 
the  amount  of  buildup  is  small  relative  to  neutron  buildup  at 
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these  ranges  because  of  the  nature  of  the  energy-dependent 
cross  sections.  Once  the  photons  have  undergone  a  signi¬ 
ficant  amount  of  scattering,  the  photon  energy  has  dropped  to 
the  point  where  the  1/E*  photoelectric  cross  section  domi¬ 
nates,  and  the  photons  are  absorbed.  Attenuation  soon  domi¬ 
nates  over  buildup  at  greater  ranges.  Attenuation  in  the 
case  of  the  photon  radiation  includes  attenuation  from 
scattering  as  well  as  from  absorption. 

2)  The  magnitude  of  the  secondary  photon  fluence  shown 
in  figure  7  is  approximately  two  orders  of  magnitude  below 
the  magnitude  of  the  primary  photon  fluence.  This  demon¬ 
strates  that  the  production  of  a  secondary  photon  is  rare 
compared  to  the  primary  photon  source  generation  rate,  which 
is  equal  to  the  neutron  source  generation  rate.  The  result 
is  a  clear  dominance  of  the  primary  photon  radiation  at  all 
ranges. 

3)  The  secondary  photons  are  generated  as  a  result  of 
neutron  interactions  with  primarily  nitrogen  and  oxygen  in 
the  atmosphere.  It  therefore  follows  that  the  secondary 
photons  should  be  coupled  in  some  way  with  the  behavior  of 
the  neutron  radiation.  In  fact,  the  secondary  photon  radi¬ 
ation  behavior  parallels  the  neutron  radiation  behavior 
because  the  neutron  radiation  acts  like  a  source  for  the 
secondary  photon  radiation.  Even  though  fluence  data  for 
ranges  greater  than  250  kilometers  was  not  calculated,  it  can 
be  inferred  from  the  results  shown  in  figure  7  that  at  very 
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long  distances  (>  2S0  kilometers),  the  secondary  photon 
radiation  would  reach  secular  equilibrium  with  the  neutron 
radiation.  The  reason  Tor  this  is  that  the  primary  photon 
radiation  at  these  distances  would  undergo  extreme  atten¬ 
uation,  and  since  the  primary  photon  radiation  is  not  self 
generating,  would  no  longer  dominate  over  the  secondary 
radiation . 

Neutron  Fluence  vs.  Angle.  Figure  9  is  a  polar  plot 
showing  the  ratio  of  SMAUG  to  MCNP  neutron  fluence  for 
selected  source  elevation  angles  versus  slant  range. 

The  MCNP  results  were  computed  using  ring  detectors.  The 
MCNP  and  SMAUG  results  were  computed  for  ranges  between  3 
kilometers  out  to  a  maximum  of  18  kilometers.  Neutron 
fluence  data  beyond  18  kilometers  was  not  computed  because  of 
the  60  kilometer  altitude  limit  in  computing  the  mass- 
integral  in  the  SMAUG-II  code.  Because  of  the  relatively 
small  number  of  data  points  for  each  slant  range  versus 
source  elevation  angle  combination,  individual  plots  of 
fluence  versus  slant  range  are  not  shown.  Figure  9  provides 
a  more  meaningful  comparison  of  SMAUG  performance  relative  to 
MCNP  performance. 

The  comparison  shows  that  SMAUG  tends  to  overestimate 
the  neutron  fluence  (except  at  90  degrees)  for  large  positive 
angles.  For  source  elevation  angles  of  60  and  90  degrees, 
the  SMAUG  to  MCNP  fluence  ratio  tends  to  decrease  slightly 
with  increasing  range,  where  the  ratio  is  close  to  1.  For 
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Figure  9  Ratio  of  SMAUG  to  MCNP  Neutron  Fluence  -for  Short 
Ranges  vs.  Angle 


source  elevation  angles  of  30  degrees  and  lower,  the  ratios 
tend  to  increase  with  increasing  range  and  decreasing 
angle.  A  maximum  ratio  o-f  1.88  is  observed  at  a  source 
elevation  angle  of  *-90  degrees  and  a  slant  range  of  18  kilo¬ 
meters.  The  -90  degree  angle  and  18  kilometer  range  combin¬ 
ation  represents  the  maximum  mass-integral  obtainable  in  this 
plot,  74. 8  g/cm* .  Figure  9  shows  a  definite  correlation 
between  the  fluence  ratios  and  increasing  density  (with 


71 


decreasin 


angle)  and  increasing  mass  range.  Since  in  every 


case  (except  at  90  degrees)  SMAUG  overpredicts  the  neutron 
fluence,  there  must  be  one  or  more  two-dimensional  effects 
which  reduces  the  neutron  contribution  to  the  free-field 
fluence  and  is  not  accounted  for  by  SMAUG-II.  This  suggests 
that  the  mechanism  of  particle  loss  in  either  the  upper  or 
lower  altitudes  is  not  accounted  for  in  SflAUG-II. 

Particle  loss  can  occur  in  several  ways* 

1)  Neutron  collisions  leading  to  scatter  away  from  the 
point  of  interest  without  further  contribution  to  the 
fluence.  A  specific  case  of  this  is  a  particle  travelling 
downward  being  scattered  in  the  upward  direction  and  leaving 
the  problem  without  further  collision. 

2)  Neutron  collisions  which  lead  to  scatter  out  of  the 
problem,  i.e.,  the  top  of  the  atmosphere,  without  further 
collision.  The  probability  for  this  mechanism  increases  with 
increasing  altitude  (decreasing  density). 

3)  Successive  scattering  of  neutrons  which  leads  to  the 
subsequent  loss  of  energy.  This  eventually  reduces  the  mean 
free  path  of  the  neutrons  to  a  point  where  the  neutrons 
diffuse  away  from  the  point  of  interest  and  can  no  longer 
contribute  to  the  fluence.  The  loss  mechanism  that  dominates 
depends  on  the  altitude  region  (and  hence  the  density),  the 
direction  of  travel,  and  the  initial  energy  of  the  neutrons. 

The  ratios  shown  in  figure  9  suggest  that  SMAUG's 
capability  to  predict  the  neutron  fluence  worsens  with 
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increasing  mass-integral ,  especially  at  -field  points  below 
the  source  altitude.  This  indicates  that  the  higher  air 
density  below  the  source  contributes  to  the  increased  attenu¬ 
ation  by  successive  scattering  of  the  neutrons  and  diffusion 
away  from  the  point  of  interest.  This  leads  to  a  significant 
loss  of  particles  in  this  region.  In  the  upper  altitude 
regions  and  the  largest  positive  angles,  SMAUG' s  ability  to 
predict  the  fluence  is  fairly  good  relative  to  the  MCNP 
results.  This  is  a  consequence  of  the  lower  air  density. 

The  collision  probability  for  neutrons  is  much  less  in  the 
upper  atmosphere  than  it  is  in  the  lower  atmosphere.  This 
means  that  most  particles  will  have  at  most  only  a  few  colli¬ 
sions  before  escaping  the  problem. 

Photon  Fluence  vs.  Angle.  Figure  10  is  a  polar  plot 
showing  the  ratio  of  SMAUG  to  MCNP  photon  fluence  for 
selected  source  elevation  angles  versus  slant  range. 

The  MCNP  results  were  computed  using  ring  detectors  for  the 
same  range  and  angle  combinations  as  shown  in  figure  9. 

The  comparison  shows  that  SMAUG  tends  to  underestimate 
the  photon  fluence  for  large  positive  angles  (except  at  very 
close  range).  For  source  elevation  angles  between  30  and  90 
degrees,  the  SMAUG  to  MCNP  fluence  ratio  tends  to  decrease 
with  increasing  range.  For  source  elevation  angles  below  30 
degrees,  the  ratios  generally  tend  to  increase  with  increa¬ 
sing  range  and  decreasing  angle.  A  maximum  ratio  of  1.74  is 
observed  at  a  source  elevation  angle  of  -90  degrees  and  a 
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Figure  10  Ratio  of  SMAUG  to  MCNP  Photon  Fluence  at  Short 
Ranges  vs.  Angle 

slant  range  of  18  kilometers. 

The  photon  radiation  behavior  parallels  the  behavior  of 
the  neutron  radiation  at  points  below  the  source,  but  not 
above  the  source.  If  the  mean  free  path  of  the  photons  was 
long  compared  to  that  of  the  neutrons,  fewer  photon  scatters 
would  occur  and  MCNP  would  predict  a  lower  photon  fluence  at 
the  point  of  interest  relative  to  the  neutron  fluence.  The 
raw  data  listed  i  i  Appendix  A,  corresponding  to  the  ratios 
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shown  in  figures  9  and  10,  substantiate  this  fact.  As  a 
consequence,  photon  leakage  would  increase  as  compared  to  the 
neutron  leakage,  especially  at  high  altitudes. 

In  a  comparison  with  SMAUG-II,  the  two-dimensional 
effects  accounted  for  by  MCNP  are  not  accounted  for  by  SMAUG- 
II.  SMAUG-II  fails  to  account  for  the  increased  particle 
leakage  in  the  lower  altitude  regions  because  particles  trav¬ 
elling  downward  have  a  high  probability  for  collision  and 
subsequent  scatter  upwards  and  out  of  the  problem.  This 
parallels  the  behavior  of  the  neutron  radiation.  The  rela¬ 
tive  differences  in  the  SMAUG  to  MCNP  photon  ratios  is  less 
than  the  corresponding  neutron  ratios.  This  suggests  that 
the  two-dimensional  effects  which  govern  the  behavior  of  the 
photon  radiation  are  less  pronounced  than  for  the  neutron 
radiation . 

Computational  Costs.  Previous  studies  have  demonstrated 
that  Monte  Carlo  transport  studies  are  particularly  suited 
for  cases  where  accurate  benchmark  results  are  needed. 
Computer  run  time  has  not  been  a  critical  issue  because  often 
only  a  few  data  points  are  needed.  In  terms  of  computer 
cost,  Monte  Carlo  transport  simulation  has  proven  to  be 
impractical  compared  to  other  methods  of  radiation  transport 
such  as  discrete  ordinates  when  large  amounts  of  data  are 
needed . 

Several  computer  simulations  using  just  the  variable 
density  version  of  MCNP  we  -e  performed  to  provide  an  under- 
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standing  of  the  added  computational  cost  associated  with 
using  ring  detectors  versus  using  no  ring  detectors.  A 
modified  MCNP  run  using  8C,000  source  neutrons  and  no  tallies 
or  detectors  for  the  25  cell  geometry  shown  in  figure  4  was 
performed.  Geometric  splitting/Russian  roulette  was  employed 
out  to  180  kilometers,  and  a  total  of  573,600  particles  were 
tracked.  The  run  time  was  153  minutes,  or  0.115  seconds  per 
source  neutron.  Another  modified  MCNP  run  was  performed 
using  the  same  geometry  and  number  of  source  neutrons,  but 
with  one  ring  detector  at  a  source  co-altitude  range  of  100 
kilometers.  The  same  number  of  total  particles  were  tracked. 
The  run  time  was  454  minutes,  or  0.341  seconds  per  source 
neutron.  As  a  result,  a  run  time  of  301  minutes  was  added  by 
using  just  one  ring  detector.  An  additional  modified  MCNP 
run  was  performed  using  two  ring  detectors  at  source  co¬ 
altitude  ranges  of  125  and  150  kilometers.  Again,  the  same 
total  number  of  source  particles  were  tracked.  The  run  time 
was  655  minutes,  or  0.491  seconds  per  source  neutron.  As  a 
result,  an  additional  201  minutes  was  added  to  the  run  time 
by  adding  one  more  ring  detector.  A  modified  MCNP  run  was 
performed  using  an  18-cell  geometry  with  18  ring  detectors 
for  ranges  out  to  250  kilometers.  For  100,000  source  neutrons 
and  geometric  splitting,  a  total  of  480,800  particles  were 
tracked.  The  run  time  was  1996  minutes,  or  1.2  seconds  per 
source  neutron, 

A  more  significant  basis  for  comparing  run  times  is  to 


compare  run  times  between  the  unmodified  version  of  MCNP  and 
the  unmodified  version  of  MCNP.  Comparisons  were  made  of  run 
times  from  previous  studies  using  ring  detectors,  specif¬ 
ically,  Culp  et  al.  C2>31-323.  An  unmodified  MCNP  run  using 
7500  source  neutrons  and  15  ring  detectors  for  the  181-cell 
geometry  shown  in  figure  2  C2:20a3  was  performed.  With 
geometric  splitting/Russian  roulette,  a  total  of  380,693 
particles  were  tracked.  The  run  time  was  580  minutes,  or 
4.64  seconds  per  source  neutron.  The  same  geometry  and  ring 
detector  locations  were  used  in  a  modified  MCNP  run.  As 
expected,  the  modified  version  of  MCNP  took  longer,  or  about 
840  minutes,  a  45  percent  increase  relative  to  the  unmod  ed 
version  of  MCNP.  This  was  a  consequence  of  the  run  time 
added  because  of  the  array  searches  required  in  the  variable 
density  model.  In  reality,  fewer  numbers  of  cells  would  be 
required  in  the  modified  version  of  MCNP  which  would  signifi¬ 
cantly  decrease  the  run  time. 

With  respect  to  the  time-intensive  array  searches  in  the 
variable  density  model,  the  run  times  from  the  modified 
version  of  MCNP  proved  to  be  very  competitive  when  the  number 
of  cells  was  reduced  in  the  modified  MCNP  code.  A  comparison 
test  was  performed  using  the  same  181-cell  geometry  as  shown 
in  figure  2  C2i20a3.  Here,  the  unmodified  MCNP  code  employed 
ring  detectors  placed  at  selected  slant  ranges  from  3  to  18 
kilometers  for  various  source  elevation  angles.  The  run  time 
results  from  the  modified  MCNP  code  were  bassd  on  the  same 
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geometry  using  only  46  cells  ar.d  identical  ring  detector 
locations.  The  comparison  test  showed  the  modified  version 
of  MCNP  ran  faster  than  the  unmodified  version  of  MCNP  by 
almost  a  factor  of  two.  No  direct  run  tine  comparisons  were 
made  using  the  unmodified  MCNP  code  at  radii  (ground  ranges) 
greater  than  100  km.  An  excessive  number  of  cells  would  have 
been  required  in  the  unmodified  version  of  MCNP  resulting  in 
an  extremely  long  run  time. 


VII.  Summary  and  Conclusions 


Summary 

Validation  of  MCNP.  Validation  of  the  modified  MCNP 
code  was  accomplished  in  two  phases.  The  first  phase  in¬ 
volved  validating  each  scale  height  factor  for  both  density 
and  mass-integral ,  and  then  validating  each  subroutine  or 
function  that  used  them.  The  scale,  height  factors  were  used 
in  the  appropriate  form  of  the  density  or  mass-integral 
fitting  equation  to  represent  the  density  variation  of  the 
atmosphere.  The  fitting  equations  produced  results  which 
were  within  1  percent  of  the  reference  data.  The  second 
phase  involved  validating  the  modified  MCNP  code  as  an  inte¬ 
grated  whole  based  on  a  two-part  approach.  In  the  first 
approach,  a  Monte  Carlo  simulation  from  a  previous  study  was 
used  to  validate,  on  a  rough  basis,  the  results  obtained  from 
the  modified  MCNP  code.  The  comparison  showed  a  maximum 
relative  difference  in  the  fluence  results  of  2  percent.  In 
the  second  approach,  a  simpler  geometry  using  fine-mesh 
spatial  grids  with  accurate  average  densities  in  each  scale 
height  region  was  incorporated  into  the  input  file.  Results 
from  the  unmodified  MCNP  code  using  this  input  file  were 
compared  to  results  from  the  modified  MCNP  code.  A  compar — 
ison  showed  a  maximum  relative  difference  in  the  fluence 
results  of  1  percent. 

MCNP  Results.  Mul tidimensional  Monte  Carlo  simulation 
studies  of  neutron-photon  transport  were  performed  within  a 
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variable  density  atmosphere.  The  calculations  were  performed 
using  the  newly-modified  version  of  MCNP  for  a  point  iso¬ 
tropic  neutron-photon  source  at  an  altitude  of  40  kilometers. 
A  collection  of  benchmark  results  using  ring  detectors  were 
obtained  for  various  angles  and  ranges.  Specifically, 
fluence  results  were  obtained  for  field  points  at  constant 
ranges  as  a  function  of  angle  and  for  field  points  at 
constant  altitudes  as  a  function  of  range.  These  results 
were  mapped  out  giving  fluence  distribution  values  as  a 
function  of  range  and  angle.  The  computer  code  SMAUG-II, 
which  uses  the  mass-integral  scaling  approximation ,  was  used 
to  generate  neutron  and  photon  fluence  results  at  the  same 
field  point  locations. 

Conclusions 

In  terms  of  the  effort  put  forth,  as  well  as  the  knowl¬ 
edge  and  insight  gained  in  modifying  a  production  Monte  Carlo 
transport  code  to  incorporate  a  variable  density  atmosphere, 
this  project  was  a  success.  The  limitations  imposed  by  using 
discrete  cells  to  represent  the  variation  in  air  density  has 
been  eliminated.  MCNP  performed  very  well  in  terms  of  saving 
computer  time,  user  overhead,  and  expanding  the  dimensions  of 
the  problem.  Validation  of  previous  studies  as  well  as 
extending  the  body  of  knowledge  in  this  area  were  accom- 
pl ished . 

Fluence  data  obtained  with  MCNP  behaved  in  an  anticipat¬ 
ed  manner,  offering  no  particular  surprises.  The  magnitude 
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of  the  neutron  fluence  was  shown  to  be  roughly  an  order  of 
magnitude  higher  than  the  corresponding  photon  results.  This 
indicated  that  neutrons  were  the  principal  form  of  radiation 
in  this  analysis.  This  also  implies  that  the  mean  free  path 
for  neutron  scatter  was  short  compared  to  that  for  photons. 

Along  the  source  co-altitude  band,  streaming  of  the 
neutron  and  photon  radiation  was  obeerveo,  indicating  reduced 
scatter  contributions  at  very  short  ranges  and  low  mass- 
integral.  For  both  neutron  and  photon  radiation,  buildup  at 
intermediate  ranges  from  scatter  increased  for  decreasing 
source  elevation  angle  and  increasing  range.  The  amount  of 
buildup  from  the  photon  radiation  was  less  pronounced  than 
that  from  the  neutron  radiation.  At  long  ranges,  attenuation 
of  the  neutron  radiation,  due  to  successive  downscatter  and 
diffusion  away  from  the  points  of  interest,  began  to  dominate 
over  buildup,  and  the  fluence  fell  off  faster  than  predicted 
by  spherical  divergence.  For  the  photon  data,  attenuation 
from  successive  scatter  and  absorption  at  low  photon  energies 
dominated  buildup  at  longer  ranges.  From  3  kilometers  out  to 
at  least  250  kilometers,  the  primary  photon  radiation  was  the 
principal  contributor  to  the  total  coupled  photon  free-field 
fluence.  It  was  shown  that  the  secondary  photon  radiation 
paralleled  the  neutron  radiation.  There  was  evidence  indi¬ 
cating  that  convergence  towards  secular  equilibrium  with  the 
neutron  radiation  at  very  long  distances  (>  250  kilometers) 
was  occurring.  This  served  as  evidence  that  the  neutron 
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r ;  ,3n  acted  as  a  source  -for  the  secondary  photon  radia- 
tic-  - 

Listed  below  are  several  factors  which  influenced  the 
MCNP  results! 

1)  neutrons  were  the  dominant  form  of  radiation; 

2)  scattering  was  the  dominant  mechanism  of  neutron- 
photon  transport  within  the  atmosphere; 

3)  two-dimensional  effects  such  as  buildup,  low-energy 
diffusion,  and  particle  escape  were  primary  factors  in  the 
relative  differences  between  MCNP  and  SMAUS-II  results. 

With  respect  to  the  approximations  used  in  tnis  analy¬ 
sis,  SMAUG-I I  predictions  of  neutron-photon  transport  charac¬ 
teristics  were  shown  to  be  accurate  to  within  a  few  percent 
at  points  within  a  few  kilometers  of  the  source.  SMAUG-I I 
neutron  results  differed  by  a  factor  of  14  along  the  source 
co-altitude  band.  The  photon  results  differed  by  a  factor  of 
6  at  this  same  elevation.  For  ranges  between  3  and  18  kilo¬ 
meters  versus  source  elevation  angle,  SMAUG-I I  again  showed 
an  inability  to  model  the  two-dimensiona 1  behavior  of  the 
neutron  and  photon  radiation  as  both  range  and  mass  range 
increased.  SMAUG-I I  had  the  greatest  difficulty  in 
predicting  the  fluence  in  regions  of  higher  density  and  over 
large  mass  ranges. 

eased  on  these  comparison  results,  there  were  several 
deficiencies  noted  with  SMAUG-I I.  First,  SMAUG-I I  could  not 
account  for  the  escape  of  particles  from  the  problem,  espe- 
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cially  out  the  top.  This  becomes  a  concern  at  higher  eleva¬ 
tions  where  particle  leakage  can  be  significant  because  of 
low  air  density  and  the  small  probability  for  particle  inter — 
action.  Also,  leakage  increased  as  a  result  of  particles 
travelling  downward,  scattering  upward,  and  leaving  the 
problem  without  further  collision.  Second,  as  a  direct 
consequence  of  the  reasons  discussed  above,  SMAuG-II  is 
inadequate  for  use  at  altitudes  above  the  original  20  kilo¬ 
meter  limit.  This  stems  from  the  inability  of  SMAUG-II  to 
account  for  multidimensional  transport  effects  inherent  in  a 
variable  density  atmosphere.  This  severely  limits  the  use¬ 
fulness  of  SMAUG-II  as  an  initial  transport  diagnostic  tool 
or  as  a  tool  for  comparison  against  other  radiation  transport 
methods.  However,  this  was  expected  since  mass-integral 
scaling  is  based  on  one-dimensional  discrete  ordinates 
results  for  a  homogeneous  system. 

Increasing  the  g,  ..metrical  dimensions  of  the  problem 
over  previously  used  geometries  in  past  Monte  Carlo  studies 
was  now  possible.  The  increase  in  run  time  as  a  consequence 
of  incorporating  a  variable  density  atmosphere  into  MCNP  was 
more  than  offset  by  the  time  saved  in  eliminating  unneeded 
spatial  cells  to  represent  the  variation  of  air  density  with 
altitude.  An  added  benefit  was  to  decrease  some  of  the  user 
overhead  in  preparing  the  input  files. 

In  terms  of  computational  costs,  using  ring  detectors 
proved  to  be  very  expensive.  Because  of  the  unique  nature  of 
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the  modifications,  the  use  of  cell  or  surface  tallies  became 
impractical.  Cells  were  still  required  but  were  so  large 
that  average  cell  or  surface  tallies  became  meaningless. 

Using  average  cell  or  surface  tallies  would  result  in  average 
fluence  values  that  do  not  accurately  represent  the  fluence 
values  within  a  large  cell,  since  the  variation  of  fluence 
across  the  cell  or  surface  is  large. 

In  comparison  to  run  times  obtained  using  an  unmodified 
version  of  MCNP,  run  times  from  the  modified  version  of  MCNP 
proved  to  be  extremely  favorable.  The  goal  of  reducing  the 
number  of  cells  required  and  at  the  same  time  offset  the 
added  run  time  caused  by  using  the  variable  density  model  was 
accomplished.  The  run  time  added  as  a  consequence  of  em¬ 
ploying  the  variable  density  atmospheric  model  was  more  than 
offset  by  the  time  saved  in  using  a  small  number  of  spatial 
cel  Is. 
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VIII.  Recommendations 


The  following  are  recommendations  for  further  study  in 
this  area.  First,  more  work  is  needed  in  developing  the 
modified  MCNP  code  so  that  it  is  applicable  to  a  wider 
variety  of  density-dependent  atmospheric  transport  problems. 
Undoubtedly,  there  are  density-dependent  parameters  and 
capabilities  employed  by  MCNP  that  were  not  modified  in  this 
project.  Capabilities  requiring  modification  might  include 
fission  heating  within  the  atmosphere,  radiation  energy 
deposition,  exponential  transf ormation  of  the  macroscopic 
cross  section  which  utilizes  path  length  stretching,  and  any 
density-dependent  dose  response  functions.  Second,  an  inves¬ 
tigation  of  alternate  geometries  and/or  variance  reduction 
schemes  to  increase  the  efficiency  of  specific  problems 
should  be  made.  Third,  alternative  array-search  algorithms 
might  improve  the  efficiency  of  the  modified  MCNP  code. 
Fourth,  MCNP  should  be  further  modified  to  allow  the  use  of 
more  than  one  material  within  the  atmosphere,  i.e,  silicon, 
tissue,  etc.  Finally,  modifications  to  allow  the  user  to 
select  either  the  variable  density  version  or  the  unmodified 
version  should  be  incorporated  into  a  single  code. 
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Appendix  Ai  MCNP  and  SMAUG  Neutron  and  Photon  Fluence  Rgsu 


MCNP  vs.  SMftUB  Neutron  Fluence  at  Selected  Altitudes  vs.  ~<a .u'j 

ALTITUDE  «  40  km 

Ground 

Range  MCNP  Neutron  Pluonce  SMAUG  Fluence 
<km)  <n/cm* /source  n)  < n /cm* ) /source  n) 


3 

1 .2046E-12 

0.0158 

1 .2551E-12 

6 

3.7991E-13 

0.0067 

A .083 IE— 13 

9 

1 .9915E-13 

0.0073 

2.2396E-13 

12 

1 . 2433E-1 3 

0.0092 

1  .4924E-13 

13 

B.3064E-14 

0.0084 

1 .0962E-13 

18 

5.8523E-14 

0.0089 

8.5201E-14 

25 

2.8385E-14 

0.0090 

5.3382E-14 

50 

3.B769E-15 

0.0120 

1 .6425E-14 

75 

9. 1952E-16 

0.0187 

6 . 3846E—1 5 

100 

2.7B91E-16 

0.0205 

2.7332E-15 

125 

1 .0528E-16 

0.0256 

1 .2609E-15 

ISO 

4.5912E-17 

0.0370 

6.2217E-16 

175 

2.3697E-17 

0.0624 

3.2513E-16 

200 

1 .4024E-17 

0.0674 

1  .7764E-16 

225 

B . 6885E- 1 8 

0.0B13 

1 .0022E-16 

250 

4 . 0998E-18 

0.0874 

5.7828E-17 

ALTITUDE  *  60  kw 

Ground 

Range  MCNP  Neutron  Fluence  SMAUG  Fluence 
(km)  (n /cm* /source  n>  (n/cm* ) /source  n) 


3 

3 . 9844E- 1 4 

0.0047 

3.8592E-14 

6 

3.7524E-14 

0.0036 

3.6795E-14 

9 

3.441  IE-14 

0.0035 

3.4212E-14 

12 

3 . 0B65E-1 4 

0,0034 

3 , 1255E-14 

15 

2.7066E-14 

0.0030 

2 . 8254E-1 4 

18 

2.3564E-14 

0.0038 

2.5417E— 14 

23 

1 .6532E-14 

0.0035 

1 .987 IE- 14 

50 

4 . 7007E-1 5 

0.0046 

9.6673E— 15 

75 

1 .5593E-15 

0.0056 

5.7184E-15 

100 

6 . 2507E-1 6 

0.0071 

3.71 99E-1 5 

125 

2.8793E-16 

0.0078 

2.5426E -15 

150 

1 .5032E-16 

0.0097 

1  . 7896E-1 5 

175 

8 . 5504E-1 7 

0.0162 

1 .2842E-15 

200 

5 . 0322E-1 7 

0.0120 

9.3478E-16 

225 

3. 1846E-17 

0.0153 

6.8B32E-16 

250 

2 . 0549E- 1 7 

0.0146 

5. 1 194E-16 
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MCNP  vs.  SMAUG  Photon  Fluence  at  Selected  Altitudes  vs.  Ranoe 

ALTITUDE  «  40  km 


Ground  MCNP  Photon  Fluence  SflAUG  Fluence 


Range 

(km) 

Primary 

< D/cm* /source  o) 

Secondary 
< D/cm* /source  n> 

Coupled 
(o/cm*  > 

3 

9.B204E-13 

0.003S 

7 . 9830E-15 

0.1120 

1 .021  IE-12 

6 

2.7361E-13 

0.0044 

4 . 01 90E—15 

0.0442 

2.8118E-13 

9 

1 . 3289E-13 

0.0063 

3. 1097E-15 

0.0589 

1 .3710E-13 

12 

7 . 9135E-14 

0.0020 

2 . 3603E-15 

0.0376 

8 . 3808E-14 

15 

5 . 2243E- 14 

0.0075 

1 . 8363E-15 

0.0355 

5 . 7683E-14 

18 

3.671BE-14 

0.0061 

1 .5446E-15 

0.0356 

4.2646E-14 

25 

1 . S368E-14 

0.0052 

8.4712E-16 

0.0541 

2 . 4706E-14 

50 

2.8510E-15 

0.0064 

2.6109E-16 

0.0596 

6.9S69E-15 

75 

7.3869E-16 

0.0095 

8 . 4884E— 17 

0.0464 

2.7902E-15 

100 

2.5194E-16 

0.0175 

3 . 3049E— 1 7 

0.0395 

1 .2695E-15 

125 

1 .0119E-16 

0.0154 

1 .4741E-17 

0.0452 

6 . 2510E-16 

150 

4.8325E-17 

0.0209 

7 . 2322E-18 

0.0508 

3.2651E-16 

175 

2-5252E-17 

0.0243 

4. 10U3E-18 

0.0783 

1 .7897E-16 

2  00 

1 .3196E-17 

0.0281 

2 . 2770E— 18 

0.1204 

1 .0218E-16 

225 

7.8418E-18 

0.0355 

1 .4547E-18 

0.1221 

6.0406E-17 

250 

S.0248E-18 

0.0459 

9.0242E-19 

0.0987 

3.6786E-17 

Ground 

Range 

(km) 

ALTITUDE  *=  60  km 

MCNP  Photon  Fluence 

Primary  Secondary 

( o/cm* /source  d)  ( o/cm* /source  n) 

SMAUG  Fluence 
Coupled 
( D/cm*  ) 

3 

2 . 7234E-14 

0.0048 

6 . 7325E— 16 

0.0243 

2.5517E-14 

6 

2.5631E-14 

0 . 0034 

6 . 5533E-16 

0.0242 

2 . 4 105E-14 

9 

2.3623E-14 

0.0056 

6 . 40B6E— 16 

0.0263 

2.2094E-14 

12 

2. 102QE-14 

0.0035 

6. 1205E-16 

0.0267 

1 . 9822E-14 

15 

1 .8506E-14 

0.0034 

5.5776E— 16 

0.0242 

1 .7553E-14 

18 

1 .6147E-14 

0.0037 

5.3139E-16 

0.0236 

1 .5447E-14 

25 

1  . 1556E-14 

0.0040 

4 . 2582E-16 

0.0222 

1 . 1461E-14 

50 

3.7173E-15 

0.0050 

1 .9149E-16 

0.0222 

4 . 8070E-15 

75 

1 .4193E-15 

0.0051 

8 . 6706E-1 7 

0.0260 

2 . 6133E-15 

100 

6 . 2569E-16 

0.0068 

4 . 2404E—17 

0.0251 

1 .6254E-15 

125 

3. 1239E-16 

0.0077 

2.229BE-17 

0.0264 

1 .0877E-15 

150 

1 .7068E-16 

0.0093 

1 .3165E-17 

0.0286 

7.6103E-16 

175 

1 .00B1E-16 

0.0119 

8.4487E-18 

0.0376 

5.4851E-16 

200 

6.4068E-17 

0.0347 

5. 9266E-1B 

0.0788 

4.0390E-16 

225 

3 . 8443E- 1 7 

0.0130 

3.7651E-18 

0.0373 

3.0235E-16 

250 

2.5629E-17 

0.0153 

2.6276E-18 

0,0393 

2.2935E-16 

MCNP  Neutron  Fluence  at  Selected  Angles  vs.  Range 


Slant 

Eiev  Angle 

Elev  Angle 

Range 

90  Degrees 

60  Degrees 

(km) 

(n/em* /source  n> 

(n/cm* /source  n) 

3 

1 . 1205E-12 

0.0072 

1 . 1359E-12 

0.0049 

6 

3 . 3085E-13 

0.0089 

3.3789E-13 

0.0049 

9 

1 .6434E-13 

0.0067 

1 .6653E-13 

0.0045 

12 

9.7B11E-14 

0.0053 

1 .00G3E-13 

0.0047 

15 

6.6532E-14 

0.0069 

6 . 7BB0E-14 

0.0047 

18 

4.8076E-14 

0.0041 

4 . 901 IE-14 

0.0052 

25 

2.70B6E-14 

0.0048 

2.7210E-14 

0.0043 

50 

7.7299E-15 

0 . 0030 

7 . 4427E-15 

0.0035 

75 

3.6B99E-15 

0.0029 

3 . 4706E-15 

0.0036 

100 

2. 1516E-15 

0.0029 

1  . 9980E-15 

0.0036 

125 

1 .4082E-15 

0.0029 

1 .2944E-15 

0.0037 

150 

9.8892E-16 

0.0029 

9.0466E-16 

0.0037 

175 

7.3331E-16 

0.0030 

6 . 6789E- 16 

0.0037 

200 

S.6126E-16 

0.0030 

5. 1379E-16 

0.0037 

225 

4.4835E-16 

0.0028 

4.0667E-16 

0.0037 

250 

3.6525E-16 

0.0028 

3.2914E-16 

0.0038 

Slant 

Elev  Angle 

Elev  Angle 

Range 

30  Degrees 

0  Degrees 

(km) 

< n/cm* /source  n) 

(n/cm* /source  n) 

3 

1 . 16B0E-12 

0.0065 

1 . 2046E-12 

0.0158 

6 

3.5122E-13 

0.0053 

3.7991E-13 

0.0067 

9 

1 .7755E-13 

0.0064 

1 .9915E-13 

0.0073 

12 

1 .0701E-13 

0.0054 

1 .2433E-13 

0.0092 

15 

7. 1744E-14 

0.0054 

8.3064E-14 

0.0084 

18 

5. 1238E-14 

0.0054 

5 . 8523E-14 

0.0089 

25 

2 . 7346E-14 

0.0049 

2 . B3B5E-14 

0.0090 

50 

6. 1621E-15 

0.0038 

3.8769E-15 

0.0120 

75 

2 . 5382E-15 

0.0036 

9. 1952E-16 

0.0187 

100 

1 .3742E-15 

0.0035 

2 . 7891E-16 

C. 0205 

125 

8.5910E-16 

0.0035 

1 .052BE-16 

0.0256 

150 

5.8759E-16 

0.0034 

4.5912E-17 

0.0370 

175 

4 . 2690E-16 

0.0034 

2.3697E-17 

0.0624 

200 

3.241 8E-16 

0.0034 

1 .4024E-17 

0.0674 

225 

2 . 5464E-16 

0.0034 

8.6885E-18 

0.0813 

250 

2 . 0509E- 1 6 

0.0034 

4 . 0998E-18 

0.0874 

eg 


MCNP  Photon  Fluence  at  Selected  Angles  vi 


Slant 

Elevation  Angle 

=  90  Degrees 

Range 

Primary  Photon 

Secondary  Photon 

<km> 

< D/cm* /source  o > 

Cr/cm* /source  n) 

3 

9.6328E-13 

0 .0058 

5.4849E-15 

0.0817 

6 

2.6051E-13 

0.0072 

2.8701E-15 

0.0607 

9 

1 . 2179E-13 

0.0056 

1.8464E-15 

0.0609 

12 

7. 1 116E-14 

0.0062 

1 . 3756E-15 

0.0486 

15 

4.6742E-14 

0.0040 

9.8375E-16 

0.0316 

18 

3.46B9E-14 

0.0404 

7 . 7951E-16 

0.0300 

25 

1.8246E-14 

0.0039 

4 . 8104E-16 

0.0046 

50 

5.0476E-1S 

0.0048 

1 .9205E-16 

0.0257 

75 

2.3378E-15 

0.0045 

1 .0260E-16 

0.0219 

100 

1 .341  IE-15 

0.0029 

6.4421E-17 

0.0222 

125 

8.6577E-1 6 

0.0029 

4. 3465E-17 

0.0220 

150 

6.0&44E-16 

0.0027 

3 . 1421E-17 

0.0234 

175 

4.4703E-16 

0.0027 

2.3698E-17 

0.0226 

200 

3.4399E-16 

0.0029 

1 .8624E-17 

0.0219 

225 

2. 7302E-16 

0.0032 

1 .5330E-17 

0.0204 

250 

2.2150E-16 

0 . 0028 

1 .0637E-17 

0.0204 

Slant 

Elevation  Angle 

-  60  Degrees 

Range 

Primary  Photon 

Secondary  Photon 

(km) 

( D/cm* /source  d> 

< o/cm* /source  n) 

3 

9.6614E-13 

0 . 0035 

6. 4555E-15 

0.0854 

6 

2.601  IE-13 

0 . 0036 

2 . 9231E-15 

0.0472 

9 

1 .2242E-13 

0.0031 

1 .89B9E-15 

0.0440 

12 

7.2109E-14 

0.0046 

1 .4052E-15 

0.0344 

15 

4 . 7538E-14 

0 . 0046 

1 .0611E-15 

0.0306 

18 

3.3B07E-14 

0.0036 

8.5263E-16 

0.0391 

25 

1 .8573E-14 

0.0030 

5 . 3234E- 1 6 

0.0038 

50 

5. 1 104E-15 

0.0024 

2. 0247E-16 

0.0227 

75 

2.3783E-15 

0.0023 

1 .0272E-16 

0.0206 

ICO 

1 .3708E-15 

0.0022 

6. 1986E-17 

0.0206 

125 

8.8B59E-16 

0.0020 

4.2257E-17 

0 . 0205 

150 

6 . 2370E- 1 6 

0.0021 

2.9772E-17 

0 . 0204 

175 

4.6208E-16 

0.0024 

2. 2368E-17 

0.0209 

200 

3.5506E-16 

0.0022 

1 .7544E-17 

0.0211 

225 

2.8173E-16 

0.0023 

1 .3799E-17 

0.0202 

250 

2 . 204  IE- 1 6 

0.0023 

1 . 1369E-17 

0.0277 
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Slant 

Elevation  Angle 

*  30  Degrees 

Range 

Primary  Photon 

Secondary  Photon 

(km> 

(o/cm1 /source  d> 

(d/cid1  /source  n) 

3 

9.6779E-13 

0.0032 

7.2728E-15 

0.0862 

6 

2,6482E-13 

0.0040 

3.2906E-15 

0.0454 

9 

1.2570E-13 

0.0070 

2.2524E-15 

0.0372 

12 

7 . 3562E-14 

0.0038 

1 . 66S8E-15 

0.0330 

15 

4 . 3604E-14 

0.0041 

1 .3008E-15 

0.0323 

ie 

3.4661E-14 

0.0044 

1 .0721E-15 

0.0349 

25 

1 . 8639E-14 

0.0029 

6.531BE-16 

0.0336 

50 

4.6263E-15 

C .0026 

2.1450E-16 

0.0276 

75 

2.0207E-15 

0.0028 

9.881  IE- 17 

0.0250 

100 

1 . 1200E-15 

0.0026 

5.5391E-17 

0.0215 

125 

7. 1297E-16 

0.0028 

3.5090E-17 

0.0210 

150 

4 . 8981E-16 

0.0024 

2 . 4437E-1 7 

0.0209 

175 

3.5764E-16 

0.0024 

1 .7813E-17 

0.0213 

200 

2.7349E-16 

0.0025 

1 .3491E-17 

0.0209 

225 

2. 1497E-16 

0.0024 

1 .0655E-17 

0.0214 

250 

1 .7355E-16 

0.0024 

8 . 5S05E- 1 8 

0.0220 

Slant 

Elevation  Angle 

■>  O  Degrees 

Range 

Primary  Photon 

Secondary  Photon 

(km) 

( D/cm* /source  d> 

(o/cm* /source  n) 

3 

9. 8204E-13 

0.0038 

7.9830E-J  5 

0.0954 

6 

2. 7361E-13 

0.0044 

4.0190E-15 

0.0442 

9 

1 .3289E-13 

0.0063 

3. 1097E-15 

0.0589 

12 

7.9135E-14 

0.0020 

2.3603E-15 

0.0376 

15 

5.2243E-14 

0.0075 

1 . 8363E- 1 5 

0.0355 

18 

3.6719E-14 

0.0061 

1.5446E-15 

0.0356 

25 

1 .B368E-14 

0 . 0052 

B.4379E-16 

0.0586 

50 

2.8510E-15 

0.0064 

2.6109E-16 

0.0596 

75 

7 . 3869E-1 6 

0.0095 

8.4884E-1 7 

0.0464 

100 

2.5194E-16 

0.0175 

3 . 3049E- 1 7 

0.0395 

125 

1 .01 19E-16 

0.0154 

1 .4741E-17 

0.0452 

150 

4 . 8325E-1 7 

0.0209 

7 . 2322E-18 

0.0508 

175 

2 . 5252E-1 7 

0.0242 

4 . 10B3E-18 

0.0783 

200 

1 .3196E-17 

0.0281 

2.2770E-18 

0.0877 

225 

7.841BE-1B 

0.0355 

1.4547E-18 

0.0961 

250 

5.0248E-18 

0.0459 

9.0242E-19 

0.0987 

90 


MCNP  vs.  SMAUG  Neutron  Fluence  at  Short 


Slant 

Range 

ikffli 

3 


6 


9 


12 


IS 


Elev 

Angle  MCNP  Neutron  Fluence  SMAUG  Fluence 
(Deo)  (n /cm* /source  n)  <n/cm« ) /source  n) 


90 

1.1205E-12 

0.0072 

1 .1869E-12 

60 

1 . 1359E-12 

0.0049 

1 . 1948E-12 

30 

1 . 1680E-12 

0.0063 

1 .2181E-12 

O 

1 .  2/'*46E-12 

0.0158 

1 . 2551E-12 

-30 

1 . 2321E-12 

0 . 0060 

1 .298BE-12 

-60 

1 .2445E-12 

0.0073 

1 .3360E-12 

-90 

1 .2375E-12 

0.0146 

1 .3509E-12 

90 

3.3085E-13 

0.0089 

3.4573E-13 

60 

3.3789E-13 

0.0049 

3.3210E-13 

30 

3.5122E-13 

0 . 0053 

3.7230E-13 

O 

3.7991E-13 

0.0067 

4.0B31E-13 

-30 

4. 1667E-13 

0.0082 

4.5777E-13 

-60 

4.S168E-13 

0.0101 

5.0572E-13 

-90 

4.6S75E-13 

0.0253 

5.2637E-13 

90 

1 .6434E-13 

0.0067 

1 .6859E-13 

60 

1 .6653E-13 

0.0045 

1 .7356E-13 

30 

1 .77S5E-13 

0.0064 

1 . 9043E-13 

0 

1 . 9915E-13 

0.0073 

2.2396E-13 

-30 

2.337SE-13 

0.0310 

2.7609E-13 

-60 

2 . 7380E-1 3 

0.0132 

3.3074E-13 

-90 

2.9770E-13 

0.0533 

3.5437E-13 

90 

9.781  IE-14 

0.0053 

1 .0064E-13 

60 

1 . 0083E-13 

0.0047 

1 .0435E-13 

30 

1 .0701E-13 

0 . 0054 

1 . 1852E-13 

0 

1 .2433E-13 

0.0092 

1 .4924E-13 

-30 

1.557BE-13 

0.0205 

2.0076E-13 

-60 

1 .8447E-13 

0.0182 

2.5032E-13 

-90 

1 .7651E-13 

0.0530 

2 . 6681 E- 13 

90 

6 . 6532E- 1 4 

0 . 0069 

6.6997E-14 

60 

6.7880E-14 

0.0047 

7.0123E-14 

30 

7. 1744E-14 

0.0054 

8. 1756E-14 

0 

8.3064E-14 

0.00B4 

1 .0962E-13 

-30 

1 .0518E-13 

0.0130 

1 . 5645E-13 

-60 

1 . 1697E-13 

0.0228 

1 .8246E-13 

-90 

1 .0646E-13 

0.0851 

1.7878E-13 

90 

4 . 8076E- 14 

0.0041 

4.7782E-14 

60 

4 .901  IE-14 

0.0052 

5.0328E-14 

30 

S. 1238E-14 

0.0054 

6.01 19E-14 

0 

5.B323E-14 

0.0089 

8 . 5201E-14 

-30 

7.3389E-14 

0.0214 

1 .2276E-13 

-60 

6 . 7539E-14 

0.0201 

1  . 1 141E-13 

-90 

4 . 6572E-14 

0.0596 

8 . 7739E-14 
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MCNP  vs.  SMAUG  Photon  Fluene  at  Short  Ranges  vs.  Ancle 


Slant 

Range 

(km) 

Elev  MCNP  Photon  Fluerce  SMAUG  Fluence 
Angle  Primary  Secondary  Coupled 
(Deo)  ( o/cm* /aource  n)  ( D/cm* /source  n)  (o/cm*) 

3 

90 

9 . 6328E-13 

0.0038 

5 . 4849E— 15 

0.0817 

1 .0040E-12 

60 

9.6614E-13 

0.0035 

6.4555E-15 

0.0854 

1.0059E-12 

30 

9.6779E-13 

0.0032 

7 . 2728E— 15 

0.0862 

1 .01 17E— 12 

0 

9.8204E-13 

0.0038 

7 . 9830E-15 

0.0954 

1 .021  IE- 12 

”30 

9.9200E-13 

0.0037 

8.3638E— 15 

0.0731 

1 . 0325E— 12 

-60 

9.9745E-13 

0 . 0050 

8.4984E-13 

0.0853 

1.0424E-12 

-90 

1 .01 12E-12 

0.0222 

8. 6136E-15 

0.0975 

1 . 0464E— 12 

6 

90 

2.6051E-13 

0.0072 

2.8701E-15 

0.0607 

2.6376E-13 

60 

2.601  IE-13 

0 . 0036 

2.9231E-15 

0.0472 

2.6S50E— 13 

30 

2.6482E-13 

0.0040 

3.2906E-15 

0.0454 

2.7107E-13 

0 

2 . 7361E-1 3 

0.0044 

4.0190E-15 

0-0442 

2.81 1SE-13 

-30 

2.8492E-13 

0.0059 

5. 1681E-15 

0.0578 

2.9525E-13 

-60 

2.9692E-13 

0.0074 

6.9139E-15 

0.0819 

3.0900E-13 

-90 

2.9705E-13 

0.0131 

6 . 8778E-15 

0.0768 

3. 1494E-13 

9 

90 

1 .2179E-13 

0.0056 

1 . 8464E-15 

0.0609 

1 .2135E-13 

60 

1 .2242E-13 

0.0031 

1 .8989E-15 

0.0440 

1 .2274E-13 

30 

1 . 2570E— 13 

0.0070 

2.2524E-15 

0.0372 

1.2751E-13 

0 

1 .32B9E-13 

0 . 0063 

3 . 1097E-15 

0.0589 

1.3710E-13 

-30 

1 .4328E-13 

0.0074 

3.9712E-15 

0.0430 

1 . 5209E-13 

-60 

1 .5293E-13 

0.0077 

5. 7999E-13 

0 . 0605 

1.6772E-13 

-90 

1  . 6045E-13 

0.0353 

5. 7562E-15 

0.0609 

1 .7443E-13 

12 

90 

7.1 1 16E-14 

0.0062 

1 .3756E-15 

0.0486 

6.9B89E-14 

60 

7.2109E-14 

0.0046 

1 .4052E-15 

0.0344 

7.0994E-14 

30 

7.3562E-14 

0.0038 

1 . 6658E-1 5 

0.0330 

7.4981E-14 

0 

7.9135E-14 

0.0070 

2.3603E-15 

0.0376 

8.3B08E-14 

-30 

8.5504E-14 

0.0067 

3.9386E-15 

0.0577 

9.8524E-14 

' 

-60 

9. 1S38E-14 

0.0134 

6. 1544E-15 

0.0698 

1 . 1240E-13 

-90 

9 . 4749E- 1 4 

0.0529 

6 . 5894E-15 

0 . 0694 

1 .1689E-13 

IS 

90 

4.6742E-14 

0.0040 

9.8375E-16 

0.0316 

4.5461E-14 

60 

4.7538E-14 

0,0046 

1 .061  IE-15 

0,0306 

4.6349E-14 

30 

4.8604E-14 

0 . 004 1 

1 .3008E-15 

0.0323 

4.96B2E-14 

0 

S . 2243E-1 4 

0.0075 

1 .8363E-15 

0.0355 

5 . 7683E- 1 4 

-30 

5.6069E-14 

0.0085 

3.2330E-15 

0.0427 

7 . 0B97E— 14 

-60 

S.4639E-14 

0.0188 

3.2842E-13 

0.0603 

7.7722E-14 

-90 

5.0426E-14 

0.0809 

7. 6478E-15 

0.0988 

7.6223E-14 

18 

90 

3 . 46B9E- 1 4 

0.0404 

7.7951E-16 

0.0300 

3. 1927E-14 

60 

3.3807E-14 

0.0036 

8.5263E-16 

0.0391 

3.2653E-14 

30 

3 . 466  IE-14 

0.0044 

1 .0721E-15 

0.0349 

3.5463E-14 

0 

3.6718E-14 

0.0061 

1 .5446E-15 

0.0356 

4.2646E-14 

-30 

3.6743E-14 

0.0095 

2 . 7201E-15 

0.0351 

5.3056E-14 

-60 

2.8615E-14 

0.0156 

4 . 1097E-15 

0.0499 

4.8631E-14 

-90 

2.0134E-14 

0.0365 

3.0505E-15 

0.0998 

4.0410E-14 

Appendix  Bt  Sample  MCNP  Input  File 


+  16N0V90  FLUENCE  (RING  DETECTORS)  AT  CO-ALTITUDE  RANGES 


101 

1 

-1 . 225E-3 

24  -25 

-23 

1MP*N=1 

IMPtP=l 

201 

1 

-1 . 225E-3 

25  -26 

-23 

IMP*N«1 

IMPlP=l 

301 

1 

-1 .225E-3 

26  -27 

-23 

!MPsN*l 

IMP«P=1 

401 

1 

-1 . 225E-3 

27  -28 

-1 

IMPiN«l 

IMPjP*1 

402 

1 

-1 . 225E-3 

27  -28 

l  -2 

IMPiN-l .21 

IMPlP-1 

403 

1 

-1 .225E-3 

27  -28 

2  -3 

lMPiN»i.44 

IMPiP*l 

404 

1 

-1 . 225E-3 

27  -28 

3  -4 

IMPiN=2.4 

IMPiPsl 

405 

1 

-1 .225E-3 

27  -28 

4  -5 

IMPtN*4.08 

IMP  i  P*  1 

406 

1 

-1 .225E-3 

27  -28 

5  -6 

IMPlN*6.75 

IMP»P«1 

407 

1 

-1 . 225E-3 

27  -28 

6  -7 

IMP*N=10.3 

IMP»P«2 

406 

1 

-1 .225E-3 

27  -28 

7  -8 

IMP*N“15.54 

IMPiP-2 

409 

1 

-1 . 225E-3 

27  -28 

8  -9 

IMPiN*20.6 

IMP»P*3 

410 

1 

-1 .225E-3 

27  -28 

9  -10 

IMP » N=29. 85 

IMP  jP=3 

411 

1 

-1.225E-3 

27  -28 

10  -11 

IMP:N=40.6 

IMP  »P«3 

412 

1 

-1 . 225E-3 

27  -28 

11  -12 

IMPlN*40 

IMP iP*5 

413 

1 

-1 . 225E-3 

27  -28 

12  -13 

I MP i N*58 . 8 

IMP  tPa5 

414 

1 

-1 .225E-3 

27  -28 

13  -14 

IMP*Ns*B9.6 

IMP  »P**8 

415 

1 

-1 .225E-3 

27  -28 

14  -15 

IMPtN«318. 1 

IMP |P=8 

416 

1 

-1 . 225E-3 

27  -28 

15  -16 

IMP»N=340.8 

IMP |P=B 

417 

1 

-1 .225E-3 

27  -28 

16  -17 

IMP»N»449 

IMP i P* 1 0 

418 

1 

-1 .225E-3 

27  -28 

17  -18 

IMPiN=637 

IMPiP«=15 

419 

1 

-1 . 225E-3 

27  -28 

18  -19 

IMP  *N* 1000 

IMP :P*40 

420 

1 

-1 . 225E-3 

27  -28 

19  -20 

IMPiN-1 

IMP iP=l 

421 

1 

-1 .225E-3 

27  -28 

20  -21 

IMP  »N*1 

IMP iP*l 

422 

1 

-1 . 225E-3 

27  -28 

21  -22 

IMP lN*l 

IMP  :P*=1 

423 

1 

-1 .225E-3 

27  -28 

22  -23 

IMP*N=1 

IMP:P«1 

501 

1 

-1 . 225E-3 

28  -29 

-23 

IMP iN=l 

IMP  »P*1 

601 

1 

-1 . 225E-3 

29  -30 

-23 

IMP*N=1 

IMP iP=l 

999 

0 

23 1-24 i 30 

IMPiN-0 

IMP*P=0 

C  23  CONCENTRIC  CYLINDERS  7  PLANES  PERPEND.  TO  Z-AXIS 


1 

CZ 

5E+5 

2 

CZ 

10E+5 

3 

CZ 

20E+5 

4 

CZ 

30E+5 

5 

CZ 

40E+5 

6 

CZ 

50E+5 

7 

CZ 

60E+5 

8 

CZ 

70E+5 

9 

CZ 

80E+5 

10 

CZ 

90E+5 

11 

CZ 

100E+5 

12 

CZ 

120E+5 

13 

CZ 

140E+5 

14 

CZ 

160E+5 

15 

CZ 

1B0E+5 

16 

CZ 

200E+5 

17 

CZ 

220E+5 

16  CZ  240E+5 

19  CZ  260E+5 

20  CZ  2S0E+5 

21  CZ  300E+5 

22  CZ  320E+5 

23  CZ  340E+5 

24  PZ  0 

25  PZ  15E+5 

26  PZ  25E+5 

27  PZ  35E+5 

28  PZ  45E+5 

29  PZ  60E+5 

30  PZ  200E+5 

MODE  N  P 

C  ISOTROPIC  NEUTRON  SOURCE  AT  (0,0,40  k licks) 

SDEF  ERG=D1  P0S*0  0  40E+5  CEL-401 
SCI  SOURCE  SPECTRUM  (GENERIC) 

SI1  H  4.14E-7  1.1 254E-6  3.059E-6  1.0677E-5  2.9023E-5  1.013E-4 
5. 0295E-4  1.2341E-3  3.3546E-3  1.0333E-2  2.1B75E-2 
2 . 4788E-2  5.2475E-2  0.1111  0.1576  0.5502  1.10B  1.827 
2.307  2.385  3.012  4.066  4.724  4.966  6.376  7.408  8.187 
9.048  10.00  11.05  12.21  12.02  13.84  14.19 
SP1  D  O  0  O  0  O  2.0226E-3  2.3974E-2 

4 . 2300E-2  7 .9823E-2  1.1337E-1  8.4647E-2  1.4145E-2 
8. 1B05E-2  7.0979E-2  3.3869E-2  9.85B4E-2  8.4961E-2 
6. 2051E-2  2.5916E-2  3.68B2E-3  2.2402E-2  2.6063E-2 
1.291BE-2  4 . 0545E-3  1.8073E-2  8.6953E-3  5.8951E-3 
6.1457E-3  7 . 8970E-3  9.4915E-3  1.6382E-2  1.736BE-2 
3 . 3667E-2  9.3298E-3 
C  MATERIAL  SPECIFICATION 

Ml  7014. 04C  -0.78  *  CONTINUOUS  NITROGEN  XSEC 

8016. 04C  -0.21  *  CONTINUOUS  OXYGEN  XSEC 

18000. 01 C  -0.01  *  CONTINUOUS  ARGON  CROSS  SECTION 

PHYSxN  14.2  *  X-SECTIONS  ABOVE  14.2  MEV  EXPUNGED 

C  TALLY  CARDS 

C  NEUTRON/PHOTON  FLUENCE  (RING  DETECTORS)  AT  250  KM 

F5Z iN  40E+5  250E+5  0 

FISZiP  40E+5  250E+5  0 

C 

NPS  100000  *  RUN  100000  HISTORIES. 

PRINT  *  PRINT  ALL  OUTPUT  FOR  EASIER  DEBUGGING 


Appendix  Ct  Sample  SMAUG-II  Input/Output  File  Listing 


Input /Output  File  Listing 


PAGE  1  SMAUG-I  I 


1 

2 

3 


( ti tle=Slant  Range  *  3  km  at  Source  Co-altitude 


) 


(yield  *1 
(read  neutron 


) 


File  * pec if ication  isi  TNS.NEV 


(0.0 

1.7368-2  1.6382-2 

(7.8970-3 

4.0545-3  1.2918-2 

(2.2402-2 

9.8584-2  3.3869-2 

(8. 1805-2 

4.2300-2  2.3974-2 

(0.0 

RSIC-36  (page  71) 

4  (norm  n«l 


1.1315-4 
9.4915-3 
6. 1457-3 
2.6063-2 
3.6882-3 
7.0979-2 
1.4145-2 
2.0226-3 
0.0 

TNS .  Neut. 

) 


5  (read  gamma 

File  specif ication  isi 


) 


9.8489-3 

) 

5.8951-3 

) 

2.5916-2 

> 

8.4647-2 

) 

0.0 

) 


EXP. GAM 


9.3298-3 
8.6953-3 
6.2051-2 
1 . 1337-1 


0.0 


(1.6496-5  1.3403-4 

8.1905-3  2.4605-2  2.7044-2 

(4.6875-2  8.1246-2 

1.0935-1  1.2896-1  4.7940-2 

(3.0055-2  2.5815-2 

1 . l*EXP(-i . l*En) .  lONovBl . 

6  (norm  g=l 

) 

7  (hob«40000 

) 

8  (alt*400C'0 

> 

9  (ran=3000 


) 

10  (go 

) 


3.0209-4 

) 

1.4082-1 

> 

1.5833-2 

) 


9.0754-4 

1.3014-1 

1.0701-2 


of  40  km 


3.3667-2 
1 .8073-2 
8.4961-2 
7.9823-2 
>.0 


2.7264-3 
1 .4655-1 
1 .0820-2 
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PAGE 


2 


S  M  A  U  G 


i 


-  I  I 


Titlei  Slant  Range  *  3  km  at  Source  Co-altitude  of  40  km 


Detonation  yield  .  .  . 

e 

1.000 

Terajoules. 

< 

.2381 

Kilotons . ) 

Detonation  altitude 

• 

4.0000E+04 

Meters . 

Prompt  neutron  yield  . 

• 

1 .000 

Neutrons . 

Primary  gamma  yield  . 

• 

1 .000 

Gammas. 

Receiver  altitude  .  . 

e 

4 . 0000E+O4 

Meters . 

Receiver  ground  range 

e 

3000. 

Meters . 

Receiver  slant  range  . 

• 

3000. 

Meters . 

Receiver  air  depth  .  . 

• 

1  .  199 

Grams/cm2 . 

A/G  Cf :  OFF,  H-A  Cfi  OFF 

,  4-Pi!  OFF. 

Neutron  mean  energy  . 

• 

1.534 

MeV. 

Neutron  fluence  .  .  . 

e 

1 .2551E-12 

N/cm2 . 

Neutron  energy  fluence 

• 

1 .9253E-12 

MeV/CM2. 

Neutron  tissue  dose 

e 

1 . 6992E-21 

Rads ( T i ) . 

Neutron  silicon  dose  . 

e 

1  . 6745E-22 

Rads(Si ) . 

Neutron  ioniz  Si  dose 

• 

1 .2777E-22 

Rads(Si)  Ion. 

1-MeV  Si-equiv  N  -fluence 

7 . 1612E-13 

N/cm2  Equiv. 

Gamma  mean  energy  .  . 

e 

.8586 

MeV. 

Gamma  -fluence  .  .  .  . 

a 

1 .021  IE-12 

G/cm2. 

Gamma  energy  -fluence  . 

a 

8.7672E-13 

MeV/cm2. 

Gamma  tissue  dose  .  . 

e 

3.9217E-22 

Rads(T i > . 

Gamma  silicon  dose  .  . 

e 

3.7927E-22 

Rads<Si ) . 

Gamma  air  dose  .... 

e 

4.1 598E-22 

Roentgens. 

N+G  Energy  -fluence  .  . 

9 

2.8020E-12 

MeV/cm2. 

N+G  Tissue  dose  .  .  . 

• 

2.0913E-21 

Rads ( Ti ) . 

N+G  Silicon  dose  .  .  . 

• 

5 . 4672E-22 

Rads ( Si ) . 

N+G  Ioniz  Silicon  dose 

• 

5 . 0704E-22 

Rads (Si)  Ion. 

N/G  Fluence  ratio  .  . 

e 

1.229 

N/G  Energy  fluence  ratio 

2.196 

N/G  Tissue  dose  ratio 

e 

4.333 

N/G  Silicon  dose  ratio 

• 

.4415 
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Explanation  of  Control  Codas: 


title® 

yield® 

read  neutron® 
norm  n* 

read  photon® 
norm  p* 


hob® 
a'i  t= 
ran* 
go 


Title  of  the  problem 

yield  of  the  source  (normally  1  Terrajoule) 
read  the  neutron  source  spectrum  file 
normalize  neutron  output  results 
(normally  1) 

read  the  photon  source  spectrum  file 

normalize  photon  output  results 

(normally  1) 

source  height  Cmetersl 

detector  altitude  Cmetersl 

detector  ground  range  Cmet*rs3 

run  problem 


Appendix  D;  Source  Code  for  INTEG  HI. BAS 


The  following  is  the  source  code  -from  INTEG_MI.BAS  which 

is  a  Microsoft  QuickBasic  V4.5  computer  code  used  to  generate 

the  mass-integral  data. 

DEFDBL  A-Z 
DEFINT  I-L,  N 

DECLARE  SUB  SIMPSON  <alt<>,  rho<),  n,  h,  MI < > ,  MI,  nbegin, 
nend) 


ver»  =  "INTEG_MI .BAS  version  1.1  of  23  AUG  90" 

* 

'  written  by!  David  L.  Monti  AFIT/GNE-91M  for  Master's  Thesis 

'  This  program  calculates  the  mass  integral  by  integrating 
a  table  of  density  values  taken  from  the  1976  U.S.  Standard 
'  Atmosphere  Tables  using  a  1/3  Simpson's  rule  integration. 

'  The  Ml's  are  then  tabulated  as  a  cumulative  distribution 
'  function  of  mass  integral  (Ml)  as  a  function  of  altitude. 


Variable  Name 


Description 


al t (  ) 
rho(  > 

MI  (  ) 
delMI (  ) 
h 
n 

nbegin 

nend 


array  of  altitude  values  Ckml 
array  of  density  values  Cg/cm^31 
array  of  integrated  mass  integral  Cg/cm'N23 
difference  between  adjacent  mass  integrals 
uniform  spacing  between  altitude  values 
number  of  data  values 

first  data  value  in  an  equispaced  group 
last  data  value  in  an  equispaced  group 


CLS 

PRINT  #2, 


a*  -  "###.#  ###.#  ###.## 

##«#. ft####**####***  ##.########" 


Open  data  file  containing  reference  atmospheric  density 
values 
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OPEN  " ATMOSREF.DAT"  ^OR  INPUT  AS  1 


OPEN  "MI . OUT”  FOR  OUTPUT  AS  2 


CALL  SIMPSQN<  al t ( ) ,  rho<>,  n,  h,  MI  <  > ,  MI,  nbegin,  nend) 


PRINT  "integrated  majs  integral  =  MI  *  100}  "Cg/cm^2D" 
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FOR  i  =  nbegin  TO  nend  STEP  2 

delMI(i)  =  ( MI ( i  +  2)  -  MI(i)>  *  100 
h  =  alt(i  +  2)  -  alt(i) 

PRINT  #2,  USING  a*;  alt(i);  altCi  +  2);  h;  rho(i);  rho(i  + 
2 > |  MI ( i >  *  100;  delMI ( i ) 

NEXT  i 

h  *  alt(62)  -  alt<61)  'every  500  meters 

nbegin  =  61  'data  every  500  meters  up  through  32  km 

nend  s  113 


'Perform  the  1/3  Simpson's  Rule  integration  on  group  3 


CALL  SIMPS0N< al t< ) ,  rho( ) ,  n,  h,  MI ( ) ,  MI,  nbegin,  nend) 

PRINT  "integrated  mass  integral  =  MI  *  100;  "Cg/cm^23“ 

FOR  i  =  nbegin  TO  nend  STEP  2 

delMI ( i )  =  ( MI ( i  +  2)  -  MI(i)>  *  100 
h  =  alt(i  +  2)  -  alt(i) 

PRINT  #2,  USING  a*;  alt<i);  alt(i  +  2);  h;  rho(i);  rho<i  + 
2);  MI ( i )  *  100;  delMI(i) 

NEXT  i 

h  =  a  1 1 ( 1 1 6 )  -  alt(llS)  'every  1000  meters 

nbegin  =115  'data  every  1000  meters  up  through  200 

km 

nend  =  201 


'Perform  the  1/3  Simpson's  Rule  integration  on  group  4 


CALL  SIMPSON  ( al  t  (  > ,  rho(  )  ,  n,  h,  MIO,  MI,  nbegin,  nend) 

PRINT  "integrated  mass  integral  *  ” ;  MI  *  100;  "Cg/cm^23" 

FOR  i  *  nbegin  TO  nend  STEP  2 

delMI ( i )  *  < MI < i  +  2)  -  MI<i>>  *  100 
h  =  al t ( i  +  2)  -  alt(i) 

PRINT  #2,  USING  a*;  alt(i);  alt(i  +  2);  h;  rho(i>;  rho(i  + 
2);  MI(i)  *  100;  delMl(i) 

NEXT  i 

h  =  alt (204)  -  alt (283)  'every  4000  meters 

nbegin  =  283  'data  every  4000  meters  up  through 

500  km 
nend  *=  356 


'Perform  the  1/3  Simpson's  Rule  integration  on  group  5 


CALL  SIMPS0N< al t< ) ,  rho( / ,  n,  h,  MI ( > ,  MI,  nbegin,  nend) 
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;  MI  *  100}  "Cg/cm^l " 


PRINT  "integrated  macs  integral  = 

FOR  i  =  n begin  TO  nend  STEP  2 

delMKi)  ■  ( MI ( i  +  2)  -  MKi))  *  100 

h  =  alt(i  ♦  2)  -  alt(i) 

PRINT  #2,  USING  a»;  alt(i)|  alt(i  +  2);  h;  rho(i>;  rhO(i 
2  > )  MKi)  *  100;  delMl(i) 

NEXT  i 

h  *  alt (359)  -  alt (358)  'every  5000  meters 

nbegin  =  358  'data  every  5000  meters  up  through 

1000  km 
nend  =  n  -  4 


'Perform  the  1/3  Simpson's  r  jle  integration  on  group  6 


CALL  SIMPSON  ( al t( > ,  rho<>,  n,  MI  (  >  ,  MI,  nbegin,  nend) 

PRINT  "integrated  mass  integral  =  MI  *  100;  "Cg/cm''23" 

FOR  i  =  nbegin  TO  nend  STEP  2 

delMI  <  i )  =  (MKi  +  2)  -  MI(i>)  *  100 
h  =  alt(i  +  2)  -  alt(i) 

PRINT  #2,  USING  at;  alt(i>;  alt(i  +  2);  h;  rho(i);  rho(i  + 
2);  MKi)  *  100;  delMKi) 

NEXT  i 
CLOSE  #2 
END 

SUB  SIMPSON  (alt<>,  rho(),  n,  h,  MK  )  ,  MI,  nbegin,  nend)  STATIC 


SUBROUTINE  SIMPSON: 

This  subroutine  performs  a  Simpson's  1/3  rule  integration 
on  a  table  of  equispacea  density  values. 

Variable  Name  Description 


rho(  ) 
a  1 1  (  ) 
n 
h 

MI 


array  of  density  values 
array  of  altitudes 
number  of  data  values 

uniform  spacing  between  altitude  values 
integrated  mass  integral 


! 


'  Integrate  using  three  data  paints  and  step  by  two 
'Use  1/3  rule 


FOR  i  »  nbegin  TC  nend  STFP  2 

MI  <  i )  «  MI  +  <h  /  3)  *  <  rho ( i )  +  4  *  rho(i  +  1)  +  rho(i  + 

2) ) 

MI  «  MI < i ) 

NEXT  i 
END  SUB 
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Appendix  Et  Source  Codes  -for  Use: — Written  Subroutines/ 

Functions 

SUBROUTINE  EQDIST (EDST ,  DTD,  ZO,  WO,  NINP) 


* 

*  THIS  SUBROUTINE  CALCULATES  2  PARAMETERS,  DEPENDING  ON 

*  INPUT: 

* 

*  1) 

t 
t 
* 

*  2) 

* 

* 

* 

* 

*  ALL  UNITS  ARE  IN  cgs,  CONSISTENT  WITH  MCNP. 

* 

* -  ■  ■■ . - . .  . - . -  ■  ■  -  ■  -  - . 

* 

*  THE  FOLLOWING  IS  A  LIST  OF  ARRAYS  AND  VARIABLES  USED  IN 

*  THIS  SUBROUTINE  AND  EXTERNAL  TO  THE  FUNCTIONS. 

* 

* 

* 

*  ZO 

*  wo 

*  EDST 

* 

*  DTD 

t  NINP 

« 

*  NR1 

* 

ItaaBaaaas  -  ■■  ■ saaaa— BsaasaanaanBaaM—— — asai  Maa—— 

include  'CM.inc' 

EXTERNAL  MAS I , DNTY , ZMAS , DM I , MI Z , EQUI VDST 
DOUBLE  PRECISION  ZO , ZD , WO ,MI CA , DENCA , DENCP , ZC 
DOUBLE  PRECISION  MASI , DNTY, ZMAS , DMI , MI Z , EQUI VDST , DTD 
DOUBLE  PRECISION  EDST,MICP,MIDA,MI INF , MIPD ,MIPDH 

C  DEFINE  MASS  INTEGRAL  AT  INFINITY  Cg/cm*  3  ( ZO  >  990  KM) 

M  J  INF  «=  1035 . 63513*  402448 

C  DETERMINE  OPTIMUM  ARRAY  SEARCH  PARAMETERS 
L*NRl/2 

IF  <  WO . GE . 0 . .AND.ZO.LT. Z I (L) ) THEN 


-  CURRENT  PARTICLE  ALTITUDE 

-  Z-DIRECTION  COSINE  RELATIVE  TO  POLAR  ANGLE 

-  EQUIVALENT  DISTANCE  TO  COLLISION  OR 

-  EQUIVALENT  MEAN  FREE  PATH  TO  BOUNDARY/DETECTOR 

-  DISTANCE  TO  BOUNDARY  OR  DETECTOR 

-  FLAG  INDICATING  WHICH  CALCULATIONS  ARE 
PERFORMED 

-  NUMBER  OF  SCALE  HEIGHT  REGIONS 


IF  THE  INPUT  VALUE  IS  THE  SAMPLED  DISTANCE  TO 
COLLISION  FROM  THE  'HSTORY.F’  SUBROUTINE,  THIS 
SUBROUTINE  RETURNS  AN  EQUIVALENT  DISTANCE  TO 
COLLISION  BASED  ON  A  VARIABLE  DENSITY  ATMOSPHERE. 

IF  THE  INPUT  VALUE  IS  THE  CALCULATED  MEAN  FREE  PATH 
FROM  THE  ' TRANSM.F*  SUBROUTINE,  THIS  SUBROUTINE 
RETURNS  AN  EQUIVALENT  MEAN  FREE  PATH  BASED  ON  A 
VARIABLE  DENSITY  ATMOSPHERE. 
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nnn  on  on  no 


Il*l 
I2-NR1 
13®  1 

ELSE I F  <  WO . GE . 0 . . AND . 20 . GE . Z I < L ) > THEN 
I1*L 
I2*NR1 
1 3*1 

ELSE I F ( WO . LT . 0 . . AND . ZO . LT . Z I <  L ) ) THEN 
1 1*L 
12*1 
1 3=- 1 

ELSEIF<WO.LT.O. . AND . ZO . GE . Z I < L >  > THEN 
1 1 *NR 1 
12*1 
1 3=- 1 

END  IF 

I F  <  N I NP . EQ . O ) THEN 

jgsa  -  - - . -s.  i— sg- .  .  =ii:.'.i— -.'Ml  ig-mgjwagjMa. 

C  PERFORM  CALCULATIONS  BASED  ON  INPUTS  FROM  * HSTORY .F ' 

C  SUBROUTINE. 


C  CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN  A 
C  HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE. 

MIPDH  *  DMI(EDST) 

C  CALCULATE  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

MICA  *  MASKZO,  II,  12,13) 

CALCULATE  REQUIRED  MASS  INTEGRAL  AT  COLLISION  ALTITUDE  IN 
A  VARIABLE  DENSITY  ATMOSPHERE. 

MI  CP  *  MIZ  <MICA,  MIPDH,  WO) 

CHECK  TO  SEE  IF  COLLISION  ALTITUDE  IS  WITHIN  PROBLEM 
LIMITS. 

IF  (MICP.GE.O. .AND.MICP.LE.MIINF)  THEN 

CALCULATE  COLLISION  ALTITUDE  GIVEN  REQUIRED  MASS  INTEGRAL 
ALONG  PARTICLE  PATH. 

ZC  *  ZMAS(MICP, II , 12, 13) 

CALCULATE  DENSITY  AT  CURRENT  PARTICLE  ALTITUDE  AND  AT 

COLLISION  ALTITUDE  FOR  CASES  WHERE  WO  IS  APPROXIMATELY 
EQUAL  TO  0. 

IF< ( ZC-ZO) .LT. 1000. )  THEN 
DENCA  *  DNTY <Z0,I1,I2,I3) 

DENCP  *  DNTY <ZC,I1,I2,J3) 

END  IF 
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C  CALCULATE  NEW  EQUIVALENT  DISTANCE  TO  COLLISION. 

EDST  =  EQUI VDST ( ZC ,  20,  WO,  MIPDH,  DENCA,  DENCP ) 
ELSE 

EDST  =  -1. 

END  IF 

ELSEIF (NINP.EQ. 1 >THEN 


PERFORM  CALCULATIONS  BASED  ON  INPUTS  FROM  'TRANSM.F' 
SUBROUTINE. 


C  CALCULATE  ALTITUDE  AT  DETECTOR  OR  BOUNDARY  CROSSING. 

ZD  «  ZO  +  WO*DTD 

IFtZD.LT.O. .OR. ZD. GT. 21 <NR2> >THEN 
WRITE ( IUO, J  < 15X , A38) ' >MSG 
EDST*-1 
RETURN 
ENDIF 

C  CALCULATE  MASS-INTEGRAL  ALONG  PARTICLE  DIRECTION  IN 
C  A  HOMOGENEOUS  ATMOSPHERE. 

C  SEA-LEVEL  ATMOSPHERE . 

MIPDH  -  DM I < DTD ) 

C  CALCULATE  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

MICA  «  MASKZO,  11,12,13) 

C  CALCULATE  MASS  INTEGRAL  AT  DETECTOR/BOUNDARY  ALTITUDE  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

MIDA  *  MASKZD, 11, 12,13) 

C  CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

IF  <ABS< ZD-ZO) .GE. 1000)  THEN 
MI PD  »  (MI DA-MICA) /WO 
ELSE 

DENCA  *  DNTY(ZO, II , 12, 13) 

DENCP  *  DNTY(ZD, 1 1 , 12, 13) 

MI PD  »  DTD*  <  DENCA+DENCP  > /2 . 

ENDIF 

C  CALCULATE  A  NEW  MEAN  FREE  PATH  BASED  ON  A  VARIABLE 
C  DENSITY  ATMOSPHERE. 

EDST  «  MIPDH  *  EDST  /  MIPD 

ENDIF 

RETURN 

END 


*nnnnnnnnnnnnn« 


FUNCTION  EQUIVDST <  ZN,  Z,  W,  MIPD,  DC A,  DCP) 


* 

C  THIS  FUNCTION  CALCULATES  A  NEW  EQUIVALENT  DISTANCE  TO 
C  COLLISION  ALONG  THE  PARTICLE  DIRECTION  WHICH  GIVES  THE 
C  SAME  HASS  INTEGRAL  AS  FOUND  IN  A  HOMOGENEOUS  SEA-LEVEL 
C  ATMOSPHERE. 

* 


VARIABLES  USED  IN  THIS  FUNCTION* 


ZN 

Z 

W 

DCA 

DCP 

MIPD 

EQUIVDST 


-  NEW  ALTITUDE  IN  A  VARIABLE  DENSITY 
ATMOSPHERE 

-  CURRENT  PARTICLE  ALTITUDE 

-  Z -DIRECT I ON  COSINE  RELATIVE  TO  THE 
POLAR  ANGLE 

-  DENSITY  AT  CURRENT  PARTICLE  ALTITUDE 

-  DENSITY  AT  NEW  ALTITUDE 

-  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION 
(HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE > 

-  EQUIVALENT  DISTANCE  TO  COLLISION  OR 
MEAN  FREE  PATH  TO  BOUNDARY/DETECTOR 


DOUBLE  PRECISION  ZN , Z ,W ,MIPD, DCA , DCP , EQUIVDST 
C 

C  IF  THE  DIFFERENCE  IN  ALTITUDES  BECOMES  SMALL  (  <  10m  >, 
C  USE  THE  AVERAGE  DENSITY  BETWEEN  THE  NEW  ALTITUDE  AND 
C  CURRENT  ALTITUDE. 

IF  (ABS(ZN-Z) .GE. 1000)  THEN 
EQUIVDST  «  (ZN-Z)/W 
ELSE 

EQUIVDST  «  MIPD  /  < O .5* ( DCA+DCP > > 

END  IF 
END 
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FUNCTION  MASI<Z, II, 12,13) 


GIVEN  AN  ALTITUDE,  FIND  THE  MASS  INTEGRAL. 

VARIABLES  USED  IN  THIS  FUNCTION: 

Z  -  GIVEN  PARTICLE  ALTITUDE 

ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 
DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 
AMI I  -  ARRAY  CONTAINING  MASS  INTEGRAL  AT  ZI 
SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL  SCALE  HEIGHT 
FACTORS 

MAS I  -  MASS  INTEGRAL  AT  GIVEN  ALTITUDE 


I  .JL  ■  ■  I  IT"  T  n  ffag3I=C=att=C=.  ■■  =^=  ■  :  tr.TS  t 

include  'CM.inc' 

DOUBLE  PRECISION  MASI.Z 
DO  10  K»I1, 12,13 

IF  (Z.GE.ZI (K> .AND.Z.LE.ZI <K+1 > )  THEN 
MAS  I  * 

Mil (KJ+DENI <K>*SHFM<K>*< 1-EXP<-(Z-ZI <K> )/SHFM(K) > ) 
END  IF 

10  CONTINUE 
END 


FUNCTION  ZMAS < MI ,11,12,13) 


GIVEN  A  MASS  INTEGRAL,  FIND  THE  ALTITUDE. 

VARIABLES  USED  IN  THIS  FUNCTION: 

MI  -  GIVEN  MASS  INTEGRAL 

SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL  SCALE  HEIGHT 
FACTORS 

ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 

AMI  I  -  ARRAY  CONTAINING  MASS  INTEGRAL  AT  ZI 

DENI  -  ARRAY  CONTAINING  DENSITV  AT  ZI 

ZMAS  -  ALTITUDE  CORRESPONDING  TO  GIVEN  MASS  INTEGRAL 


include  'CM.inc' 

DOUBLE  PRECISION  MI , ARG, ZMAS 
DO  30  KsI 1 , 1 2 , 13 

IF  ( MI . GE . AMI  I (K) . AND. MI .LE.AMII(K+1>)  THEN 
ARG  «  LOG< 1  +  <  AMI  I < K  >  -  MI )  /  (DENIOO  * 

SHFMOO  )  ) 

ZMAS  «  Z I  ( K >  -  SHFMOO  *  ARG 
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30 


) 

E^ND  IF 
CONTINUE 
END 


FUNCTION  DNTY (Z, II, 12, 13) 


* 

C  GIVEN  AN  ALTITUDE,  FIND  THE  DENSITY. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTION: 

C 

C  Z  GIVEN  PARTICLE  ALTITUDE 

C  ZI  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 

C  DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 

C  SHFD  -  ARRAY  CONTAINING  DENSITY  SCALE  HEIGHT 

C  FACTORS 

C  DNTY  -  DENSITY  CORRESPONDING  TO  GIVEN  ALTITUDE 

* 

♦  ..IT1,1:  '1U.-J1  •  1 tr  •=u-3»-g— i  -T. 1  .r-.  eawr  .a:  aa  BM—ggaa  ■■■■,  assasarn-.-igTSsaacara  i.'J.'.str. — :~«.r 

.include  'CM.inc' 

DOUBLE  PRECISION  Z,DNTY 
DO  10  K=I1  , 12,13 

IF  (Z.GE.ZKK)  .AND.Z.LE.ZI<K+1>  >  THEN 

DNTY  *  DENI (K )  *  EXP< -< Z-Z I ( K ) )  /  SHFD(K) ) 

END  IF 

10  CONTINUE 
END 


FUNCTION  DMI (DIST) 


CALCULATE  MASS  INTEGRAL  BETWEEN  TWO  POINTS  IN  A 
SEA-LEVEL  HOMOGENEOUS  ATMOSPHERE. 

VARIABLES  USED  IN  THIS  FUNCTION: 

DIST  -  DISTANCE  TO  COLLISION  OR  MEAN  FREE  PATH  TO 
BOUNDARY/DETECTOR 

DMI  -  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION 


DOUBLE  PRECISION  DIST, DM I 

CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  PATH. 
DMI  =  1 .225E-3*DIST 

END 


FUNCTION  MI Z ( MIP ,  DMI,  W) 
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FUNCTION  MIZtMIP,  DMI,  W> 


* 

C  CALCULATE  MASS  INTEGRAL  AT  NEW  ALTITUDE  IN  AN  EXPONENTIAL 
C  ATMOSPHERE. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTIONi 

C 

C  MIP 

C 

C  DMI 

C 

c  w 

c 

C  MIZ 

C 

* 


DOUBLE  PRECISION  MIP, DMI ,W, MIZ 

CALCULATE  MASS  INTEGRAL  AT  NEW  ALTITUDE. 
MIZ  =  W*DMI  +  MIP 
END 


-  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE 
IN  A  VARIABLE  DENSITY  ATMOSPHERE 

-  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN  A 
HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE 

-  Z-DIRECTION  COSINE  RELATIVE  TO  THE  POLAR 
ANGLE 

-  MASS  INTEGRAL  AT  NEW  ALTITUDE  IN  A 
VARIABLE  DENSITY  ATMOSPHERE 


BLOCK  DATA  ATINIT 


THIS  DATA  BLOCK  IS  USED  TO  SUPPORT  THE  VARIABLE 
DENSITY  ATMOSPHERIC  MODEL.  INITIALIZE  COMMON  /ATMCOM/ 
WITH  VALUES  OF  ATMOSPHERIC  PARAMETERS. 


*, 

* 

C 

c 

c 

c 

c 

c 

c 

* 


include  'CM.inc 


SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL 
FACTORS 

SHFD  -  ARRAY  CONTAINING  DENSITY  SCALE 
ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE 
REGIONS 

DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 
AMI  I  -  ARRAY  CONTAINING  MASS  INTEGRAL 


SCALE  HEIGHT 

HEIGHT  FACTORS 
OF  VERTICAL 


AT  ZI 


bb— — — eaamgeas — — — 

* 

DATA  SHFM/1034124. 022372781 , 10099S1 .924149343, 

1  986693 . 1 726 1 9466 , 964303 . 332721 5038 , 

2  941135.9145785341 ,888092.312585042, 

3  692868.8443723,637442.9797604, 
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4 

5 

6 
7 

e 

9 

+ 


+ 

+ 

+ 

+ 


+ 

DATA 

1 

2 

3 

4 

5 

6 

7 

8 
9 


+ 

+ 

+ 

+ 

+ 

DATA 

1 

2 

3 

4 

5 

DATA 

1 

2 

3 

4 

5 

6 

7 

8 
9 

DATA 

1 

2 

3 


626 1 92 . 9939628 , 639736 . 906 17414, 
47557.78050254 ,710362 . 42751 12881 , 
838211.2645959415,776503.932153373, 
683277.3511724839,612608.3454716644, 
557270 . 926 1 50725 , 568864 . 8492263636 , 
619813.9513853799,935109.6091771976, 

1 26721 7.41 5038409 , 1 578526 . 59297959 , 

1 870150 . 85627943 , 22 1 9020 . 933057442 , 

27 1 4282 . 1 34795458 , 3257224 .611 376 1 4 , 

399 1 584 . 358044663 , 4584848 . 61083246 , 

5 1 62659 . 76 1 537237 ,5791446. 057543592 , 

550 1 266 . 555948576 , 6836368 . 561027 1 77 , 
7834040 . 063075421 , 98 1 4372 . 605576273 , 

1 3917805 . 4965649 , 1 8877525 .61508398/ 

SHFD/ 1 03039 1 . 72606848B , 1 006927 . 055016296 , 

983 1 53 . 5607496358 , 960534 .1411 264988 , 
937232.106253662,866352. 1969023093, 
664086 . 7633898 1 56 , 637638 . 46529 , 

627630 . 265 119, 642604 . 7540969 , 

654589 . 3995756 , 7360 1 1 . 6057506 1 85 , 
834195.106975831 ,758287.3616637037, 
666098 . 1 02745308 , 592758 . 528 1 268486 , 
553227 . 469727463 1 , 5704 1 3 . 0947863342 , 

678 1 76 . 4068898259 , 997277 . 94 1 535646 , 

1 324262 . 26009301 3 , 1 632165 . 67546257 1 , 
1919412. 5389299 , 23 1 3392 . 496297635 , 

27974 1 2 . 470495623 , 3408285 . 2B6706647 , 

4 1 1 2545 . 885579616 , 4680457 . 5527 1 1 085 , 

53 1 1 463 . 853509 1 05 , 5890652 . 50294756 , 

6360 111.1 583 18114, 6989325 . 5647 1 1 994 , 

8 1 62740 . 323944079 ,1041 9949 . 4345468 , 
14188168. 296 14109,1 9869924 . 24084625/ 

Z 1 /0 . , 1 . D+5 , 2 . D+5 , 3 • D+5 , 4  « D+5 , 5 • D+S , 1 . D+6 , 

1 . 5D+6 , 2 . D+6 , 2 . 50+6 , 3 . D+6 , 4 . D+6 , 5 . D+6 , 

6 • D+6 , 7 . D+6 , 8 . D+6, 9 . D+6 , 1 . D+7 ,1.1 D+7 , 

1 . 2D+7 , 1 . 3D+7 , 1 . 4D+7 , 1 . 5D+7 , 1 . 6D+7 , 1 . 8D+7 , 
2 . 0D+7 , 2 . 4D+7 , 2 . 8D+7 , 3 . 2D+7 , 4 . OD+7 , 4 . 8D+7 , 
5 . 6D+7 , 6 . 4D+7 , 7 . 2D+7 , 8 . OD+7 , 8 . BD+7 , 9 . 9D+7/ 
DENI / 1 . 2250-3 ,1.111 7D-3 , 1 . 0066D-3 , 9 . 0925D-4 , 

8 . 1 935D-4 , 7 . 3643D-4 ,4.1351 D-4 , 1 . 9476D-4 , 
8.8910-5,4 .0084D-5, 1 .8410-5,3.99570-6, 

1 . 0269D-6 , 3 • 0968D-7 , 8 . 2829D-8 , 1 . 8458D— 8, 
3.416D-9,5.604D-10,9.708D-11 ,2.2220-11 , 

8. 152D-12, 3.831 D- 12,2.0760-1 2, 1.233D-12, 
5. 194D-13, 2.541 D-13, 7.8580-14,2. 971D-14, 

1 . 2640- 14,2. 803D- 15,7. 2080- 16,2. 049D- 1 6 , 

6 . 523D-1 7 , 2 . 448D- 17,1.1 36D-14 ,6. 464D-18, 
3.7160-18/ 

AMI I/O. ,116.7635,222.607167,318.334333, 

404.704533,482.4368,763.994467,91 1.2723, 
978.75926666666669,1009.379620, 

1 023 . 286286666667 , 1 032 . 662946666667 , 
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1034 . 806793333334 , 1035 . 406480000001 , 
1 035 . 580609766667 , 1 035 . 624 1 07866667 , 
1035. 633205 1 46667 , 1 035 . 634792366667 , 
1035.635056196001 , 1035.635104380667, 
1 035 . 635 1 1 8027400 ,1035. 635 1 23665300 , 
1035.635126503134,1035.635128111100, 
1035.6351 29736200 , 1 035 . 635 1 3047 1237, 
1 035 . 635 1 3 1 056504 , 1 035 . 635 1 3 1 2550 1 7 , 
1 035 . 635 1 3 1 334304 , 1 035 . 635 1 3 1 385704 , 
1035.635131397859,1035. 635 1 3 1 400898 , 
1035.635131401864,1035.635131402191 , 
1035.635131402325,1035.635131402394, 
1035.635131402448/ 
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Neutron  Spectrum  Energy  Groups 


Group 

Energy 

Band 

1  CM»VD 

1 

16.91 

to 

19.64 

2 

14.92 

to 

16.91 

3 

14.19 

to 

14.92 

4 

13.84 

to 

14.19 

5 

12.82 

to 

13.84 

6 

12.21 

to 

12.82 

7 

11.05 

to 

12.21 

8 

10.00 

to 

11.05 

9 

9.048 

to 

10.00 

10 

8.187 

to 

9.048 

11 

7.408 

to 

8.187 

12 

6.376 

to 

7.408 

13 

4.966 

to 

6.376 

14 

4.724 

to 

4.966 

15 

4.066 

to 

4.724 

16 

3.012 

to 

4.066 

17 

2.385 

to 

3.012 

1  8 

2.307 

to 

2.385 

19 

1.827 

to 

2.307 

20 

1.108 

to 

1 .827 

2' 

0.5502 

to 

1 . 108 

22 

0. 1576 

to 

0.5502 

23 

0.1111 

to 

0.1576 

24 

5.2475E-? 

to 

0.1111 

25 

2.47B8E-2 

to 

5.2475E-2 

26 

2.1B75E-2 

to 

2.4788E-2 

27 

1 . 0333E-2 

to 

2. 1875E-2 
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Neutron  Spectrum  Energy  Groups  <Cont> 


Group 

Energy  Bend 

CMeVJ 

28 

3.3546E-3 

to 

1 . 0333E-2 

29 

1.2341E-3 

to 

3.3546E-2 

30 

5.8295E-4 

to 

1 . 2341E-3 

31 

1.0130E-4 

to 

5.8295E-4 

32 

2. 9023E-5 

to 

1 .0130E-4 

33 

1 .0677E-5 

to 

2.9023E-5 

34 

3.0590E-6 

to 

1 .0677E-5 

35 

1 . 1254E-6 

to 

3.0590E-6 

36 

4. 1400E-7 

to 

1 . 1254E-6 

37 

1 .0000E-1 1 

to 

4. 1400E-7 
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Photon  Spectrum  Energy  Groups 


Group 

Energy 

Band  CMeVD  f 

1 

10.00 

to 

14.00 

2 

8.00 

to 

10.00 

3 

7.00 

to 

a> 

• 

o 

o 

4 

6.00 

to 

7.00 

S 

5-00 

to 

6.00 

6 

4.00 

to 

5.00 

7 

3.00 

to 

4.00 

8 

2.50 

to 

3.00 

9 

2.00 

to 

2.50 

10 

1  .50 

to 

2.00 

11 

1  .00 

to 

1.50 

12 

0.70 

to 

1.00 

13 

0.45 

to 

0.70 

14 

0.30 

to 

0.45 

15 

0.15 

to 

0.30 

16 

0.10 

to 

0.  15 

17 

7.00E-2 

to 

0.  10 

18 

4.50E-2 

to 

7.00E-2 

19 

3.00E-2 

to 

4.50E-2 

20 

2.00E-2 

to 

3.00E-2 

WBM 

1 .00E-2 

to 

2.00E-2 

1  14 
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